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Streszczenie 


Zrozumienie procesu powstawania tzw. ciepłego deszczu, (tzn. formowania się opadu 
atmosferycznego wskutek złożonych oddziaływań hydrodynamicznych zawiesiny kropel 
iatmosfery) jest wciąż nierozwiązanym problemem fizyki atmosfery, którego wyjaśnienie 
miałoby wielkie znaczenie dla współczesnej meteorologii i klimatologii. Trudności teoretycz- 
ne i brak wiarygodnych informacji doświadczalnych powodują konieczność stosowania 
w modelach atmosferycznych bardzo uproszczonych parametryzacji. Powszechnie przyjmuje 
się, że krople deszczowe powstają w wyniku zderzeń kropel chmurowych. Dla rozwiązania 
stochastycznego równania opisującego ewolucję rozmiarów zespołu kropel wskutek zderzeń 
konieczna jest znajomość statystyki prędkości kropel w polu turbulencji wewnątrz chmury. 
Pomiarów tego rodzaju statystyk nie przeprowadzono wcześniej, co więcej w światowej lite- 
raturze naukowej nie ma żadnych prac pomiarowych dokumentujących własności drobnoska- 
lowej turbulencji w chmurach. Celem naukowym pracy było wypełnienie tej luki i uzyskanie 
danych pomiarowych o oddziaływaniach hydrodynamicznych i termodynamicznych w struk- 
turze chmurowej poprzez analizę statystyczną pól prędkości kropel w laboratoryjnym modelu 


chmury. 


W ramach pracy skompletowany został układ eksperymentalny umożliwiający badania 
laboratoryjne procesu turbulentnego mieszania różnych mas powietrza z udziałem parujących 
kropelek wodnych. Głównym elementem układu była szklana komora o wymiarach 1,8m x 
Im x Im, w której zachodził badany proces mieszania obserwowany przez kamerę CCD, 
umieszczoną na zewnątrz komory. Zastosowanie anemometrii obrazowej (z angielskiego PIV 
— Particle Image Velocimetry) pozwoliło na bezinwazyjny pomiar pól prędkości na siatce 
o rozdzielczości poniżej skali Kołmogorowa i utworzenie bazy danych pomiarowych, która 
posłużyła później do statystycznej analizy i wyznaczenia między innymi takich wielkości jak 
funkcje struktury, korelacje, współczynnik dyssypacji lepkiej, skala Kołmogorowa. Wykona- 
nie eksperymentów w kontrolowanych warunkach umożliwiło znalezienie wpływu początko- 
wych parametrów termodynamicznych na wielkości charakteryzujące turbulencję. Porówna- 
nie danych doświadczalnych z analizą modelu mieszania z udziałem parujących kropelek 
sugeruje wpływ parowania na strukturę turbulencji wskutek zmian sił wyporu. Badania wyko- 
nane w komorze dla nieparujących kropelek z syntetycznego oleju wskazują na znaczne 
różnice struktury przepływu z parowaniem lub bez. Fakt ten jest dodatkowym argumentem 
potwierdzającym hipotezę o znaczącym wpływie parowania kropelek chmurowych na struk- 


turę przepływu w małych skalach, istotnych dla procesów koalescencji kropel, poszerzenia 


ich widma i tworzenia się opadu. Również rezultaty przeprowadzonych eksperymentów 
numerycznych oddziaływania parujących kropelek ze strukturą wirową potwierdzają tę hipo- 
tezę. 


Wykonane badania w komorze chmurowej pozwoliły na stworzenie unikalnej bazy danych 
empirycznych do walidacji modeli numerycznych procesu oddziaływań struktur turbulent- 


nych z zawiesiną kroplowa w obecności przemian fazowych. 


Effect of cloud water on small-scale turbulence - 


laboratory model 


Abstract 


Our understanding of complex hydrodynamical interactions in water droplets — air suspen- 
sion is still insufficient to explain warm rain formation problem - formation of precipitation in 
clouds without ice crystals. It is unresolved problem of atmospheric physics and progress in 
this area would give great benefits for meteorology and climatology. Difficulties in theoretical 
modelling of such complex systems and lack of reliable experimental data causes that rough 
parametrization should be applied in atmospheric models. It is assumed that rain drops arise 
as a result of water droplets collisions. Statistical properties of turbulent motion inside a cloud 
are needed to resolve stochastic equations of droplets collisions what leads to their sizes evol- 
ution. Measurements of such statistics are not available, moreover there are no reports docu- 
menting measurements of small-scale turbulence in clouds. The aim of this work was to fill in 
the gap in experimental data by providing documentation of hydrodynamical and thermody- 
namic interactions in cloud suspension based on statistical analysis of droplets velocities in 


laboratory model of cloud. 


The experimental setup was constructed to investigate process of air mixing in the pres- 
ence of evaporating water droplets in the laboratory. The main part of the setup was glass 
walled chamber 1.8m x 1 m x Im. In this chamber the process of mixing was analysed using 
CCD camera and pulsed laser illumination technique. Application of Particle Image Veloci- 
metry enabled us to measure velocity field with resolution down to Kolmogorov scale without 
disruption of investigated flow. Statistical analysis of collected data enabled evaluation of 
such characteristics of turbulence as the structure function, the turbulence dissipation rate and 
the Kolmogorov scale. Parameters describing conditions of each experiment were controlled 
what enabled to investigate influence of thermodynamic conditions on evaluated characterist- 
ics of turbulence. The experimental data and analysis of mixing model taking into account 
droplets evaporation indicate influence of buoyancy fluctuations on the turbulent flow proper- 


ties. It is confirmed by control experiments with nonevaporative droplets. 
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Cel i zakres pracy 


Ważnym problemem współczesnej fizyki chmur jest powiązanie procesów zachodzących 
w małych skalach (rzędu milimetrów i mniej) z makroskopowymi właściwościami chmur 
takimi jak ewolucja widma kropel, transfer radiacyjny, czy prędkość formowania opadu. 
Pomimo prowadzonych od kilkudziesięciu lat badań nie zostały do końca wyjaśnione zacho- 
dzące w chmurach zjawiska jak na przykład powstawanie ciepłego deszczu, czyli opadu 
atmosferycznego powstającego bez zarodników lodu. Przeprowadzone badania laboratoryjne 
mają pomóc w wyjaśnieniu roli lokalnej turbulencji przepływu na koalescencję kropel chmu- 


rowych i ewolucję ich widma. 


Celem niniejszej pracy jest analiza procesów zachodzących w małych skalach w czasie 
turbulentnego mieszania powietrza chmurowego z czystym powietrzem w warunkach labora- 
toryjnych i zbadanie wpływu parowania kropel na charakterystyki przepływu. Za powietrze 
chmurowe uważamy nasycone powietrze (wilgotność względna 100%) wypełnione kropelka- 
mi. Wykonanie eksperymentów w laboratorium pozwala na kontrolowanie istotnych parame- 
trów fizycznych procesu mieszania takich jak wilgotność czystego powietrza, temperatura, 
wielkość kropel, wodność. Pomiary pól prędkości w kontrolowanych warunkach dostarczą 
dużej liczby danych do analizy statystycznej właściwości przepływu. Pozwoli to na znalezie- 
nie wartości, które charakteryzują lokalną turbulencję w procesie mieszania, takich jak funk- 
cje struktury, dyssypację lepką, skalę Kołmogorowa czy szerokości rozkładów fluktuacji 
prędkości. Uzyskane rezultaty umożliwią weryfikację istniejących modeli procesów tworze- 
nia się opadu i będą stanowić cenną bazę doświadczalną pozwalającą na testowanie i weryfi- 
kację tworzonych w różnych ośrodkach badawczych zaawansowanych modeli numerycznych 
procesów aglomeracji kropel deszczowych. Pozwoli to na weryfikację istotnej w modelach 


parametryzacji stosowanej do przewidywania opadów. 


Dla zrealizowania sformułowanego wyżej celu został zbudowany układ eksperymentalny, 
składający się głównie z komory chmurowej o wymiarach 1m x 1m x 1,8m, układu oświetla- 
jącego, przesuwu i układu rejestrującego obrazy przepływu. W prowadzonych badaniach 
układ rejestracyjny i układ formujący wiązkę były przymocowane do ruchomej, pionowej 
prowadnicy sterowanej komputerem, co umożliwiało pomiary na różnych wysokościach w 


celu zbadania własności przepływu w różnych etapach mieszania. 
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Głównym narzędziem eksperymentalnym użytym w pracy jest anemometria obrazowa 
(z angielskiego PIV — Particle Image Velocimetry). Jest to technika pozwalająca na pomiar 
pola prędkości w przepływie, wykorzystująca obrazy poruszających się z przepływem cząstek 
wskaźnikowych (posiewu). Pomiary do opisywanych badań wykonano dla rzutów pola pręd- 
kości na dwuwymiarową płaszczyznę pomiaru wyznaczoną przez płaską wiązkę lasera oświe- 
tlającego przepływ. Metoda PIV pozwoliła na zebranie dużej liczby danych do analizy staty- 


stycznej właściwości przepływu. 


Wykonanie planowanych badań było związane z koniecznością opracowania specjalistycz- 
nego oprogramowania dostosowanego do specyficznych warunków eksperymentu i służącego 
do: 


- analizy obrazów zebranych w czasie eksperymentów i wyznaczenia pól prędkości, 
— pomiaru wielkości kropelek, 
- analizy zmierzonych pól prędkości w celu oszacowania charakterystyk turbulencji. 


Praca zawiera również analizę teoretyczną modelu mieszania z uwzględnieniem efektu 
parowania (diagram mieszania), której wyniki zostały odniesione do pomiarów w celu zbada- 
nia wpływu przemian fazowych na przepływ. Wyniki eksperymentalne uzupełniono analizą 
numeryczną procesu parowania pojedynczej kropli wody przemieszczającej się w mieszaninie 


powietrza i pary wodnej. 
W ramach przeprowadzonych w tej pracy badań: 
- wyznaczono charakterystyki mieszania powietrza chmurowego z czystym powietrzem, 


- zbudowano stanowisko pozwalające na badanie struktur przepływu powietrza chmuro- 


wego w warunkach laboratoryjnych, 


— Opracowano i wykonano nowy algorytm analizy pól prędkości, zaadaptowany do 


potrzeb eksperymentów z chmurami, 


- wyznaczono podstawowe charakterystyki struktur turbulentnych w typowych przepły- 


wach powietrza chmurowego, 
- znaleziono wpływ przemiany fazowej na charakterystyki turbulentne, 


— przeprowadzono analizę numeryczną kropli w kontekście wpływu procesu mieszania 


warstw powietrza chmurowego. 
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Ważną częścią pracy jest opracowanie specjalnego oprogramowania, służącego do pomia- 
rów pól prędkości i wielkości kropelek. Oprogramowanie to znalazło zastosowanie również 
w innych projektach badawczych [1]. Podobnie, stworzona w ramach pracy baza danych 
empirycznych procesu mieszania w komorze chmurowej jest już wykorzystywana do walida- 
cji symulacji numerycznych generacji fluktuacji turbulentnych poprzez efekty wypornościo- 


we związane z procesami parowania i kondensacji w chmurowej zawiesinie kroplowej [2]. 
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Obecność chmur w atmosferze ziemskiej odgrywa kluczową rolę w bilansie energetycz- 
nym i hydrologicznym. Stanowią więc one istotny element kształtujący oblicze naszej planety 
niezbędny w utrzymaniu odpowiednich warunków do istnienia życia. Poznanie mechanizmów 
odpowiedzialnych za cykl powstawania i skraplania się struktur chmurowych ma więc 


znaczenie nie tylko poznawcze ale też praktyczne. 


Jednym z intrygujących problemów we współczesnych badaniach nad chmurami i klima- 
tem jest powiązanie procesów zachodzących w małych skalach z makroskopowymi właści- 
wościami chmur takimi jak ich czas trwania, rozmiar, wydajność opadu i transfer radiacyjny. 
Procesy drobnoskalowe są jednak nadal słabo poznane. Wiedza w tej dziedzinie jest jeszcze 
zbyt mała by wyjaśnić takie zjawiska jak poszerzenie widma kropel i powstanie tak zwanego 


ciepłego deszczu — czyli opadu z chmury, w której nie ma kryształków lodu. 


Wzrost kropel przez kondensację nie jest wystarczająco wydajnym procesem, aby mógł 
powstać opad. Duże krople mogą tworzyć się w wyniku zderzania i koalescencji, jednak aby 
powstała kropla deszczowa o średnicy 1mm potrzeba zderzyć milion kropelek chmurowych 
o średnicy 10um. Do tego niezbędne jest występowanie fluktuacji rozkładu przestrzennego 


kropel i względnych prędkości. 


W ostatnich latach pojawiła się znaczna liczba prac pokazująca, że to właśnie turbulencja 
ma istotny wpływ na przestrzenny rozkład kropel ([3], [4], [5], [6], [7]). Tworzenie się obsza- 
rów zwiększonej koncentracji kropel i powstawanie silnych przyspieszeń związanych 
z niejednorodnością i intermittencją turbulentnego pola prędkości ułatwia proces zderzania się 
kropel co skutkuje poszerzeniem widma wielkości kropel ([8], [9], [10], [11], [12], [13]) 
i w konsekwencji powstaniem opadu ([14], [15], [16]). 


Turbulencja chmurowa powiązana z mieszaniem się dwóch mas powietrza może powodo- 
wać dodatkowe efekty które wymagają dogłębnego zbadania [17]. Woda w atmosferze, czy to 
w postaci pary, czy w postaci ciekłej zawartej w kropelkach nie jest pasywna. Szczególnie 
istotne są tutaj procesy przemian fazowych. Para wodna gromadzi pokaźne zasoby energii - 


tak zwane ciepło utajone. Jest to energia, która wyzwala się podczas kondensacji. Dla wody 
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jest to 23:10” Jkg”. Jest to niezwykle istotne dla procesów atmosferycznych. Łatwo można 
oszacować że energia uzyskana z kondensacji danej masy pary wodnej jest równa energii 
potencjalnej, gdy daną masę wyniesie się na wysokość 23-10* m. Para wodna jest dzięki temu 
„paliwem” napędzającym konwekcję atmosferyczną. Efekt tego procesu może być wyraźnie 
obserwowany na przykład w chmurach burzowych, czy cyklonach, gdzie kondensacja pary 
wodnej wyzwala energię, która powoduje powstanie potężnych prądów wstępujących wsku- 


tek zwiększenia wyporu ogrzanych mas powietrza. 


W przypadku prezentowanych tu badań przedmiotem zainteresowania był proces przeciw- 
ny, mianowicie parowanie ciekłej wody zawartej w kropelkach w procesie mieszania zawiesi- 
ny chmurowej z czystym powietrzem. Mieszające się masy powietrza tworzą skomplikowane 
struktury o różnorodnych kształtach i skalach [18]. Prowadzi to do występowania w objętości 
przepływu obszarów zróżnicowanych pod względem wielkości, kształtu, koncentracji kropel, 
wilgotności i innych parametrów termodynamicznych. Krople wpadające w obszary o małej 
wilgotności odparowują, czemu towarzyszy pobranie z otoczenia ciepła parowania. Ochło- 
dzone wskutek tego objętości powietrza charakteryzują się mniejszą wypornością. Odparowy- 
wanie wody powoduje więc zmiany parametrów termodynamicznych, które w warunkach 
grawitacyjnych pociągają zmiany siły wymuszającej związanej z wyporem. Ponieważ miesza- 
nie odbywa się w sposób chaotyczny, można również spodziewać się chaotycznego rozkładu 
obszarów zróżnicowanych pod względem sił wyporu na nie działających. Można więc 
spodziewać się, że w mieszaniu turbulentnym generowana będzie dodatkowo siła wymuszają- 


ca o nieregularnym rozkładzie w przestrzeni i czasie. 


Powstawanie dodatkowej siły wymuszającej może dodatkowo modyfikować lokalny prze- 
pływ turbulentny, a tym samym wpływać na procesy aglomeracji kropel. Przyjmuje się, że 
źródłem energii przepływu turbulentnego jest ruch średni, który wskutek działania sił inercyj- 
nych generuje duże wiry, w których gromadzi się energia turbulencji. Siły bezwładności 
powodują, że wiry te rozpadają się na mniejsze. Rozpad kolejnych wirów tworzy tak zwaną 
kaskadę energii w której energia kinetyczna przepływu przekazywana jest do coraz mniej- 
szych skal. Kaskada ta kończy się na wirach na tyle małych, że siły związane z lepkością 
płynu powodują dyssypację energii kinetycznej turbulencji. Ta kaskada wirów jest postacią 
chaosu czasowo — przestrzennego, w którym nie ma żadnej wyróżnionej struktury, a znacze- 
nia nabiera szeroki zakres skal. Postulowany przez Kołmogorowa [19] schemat spektralny 


przewiduje, że w zakresie od skali pompowania po skalę dyssypacji występuje zachowanie 
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izotropowe, w którym widmo mocy fluktuacji składowej prędkości wzdłuż dowolnej osi 


zmienia się proporcjonalnie do potęgi -5/3 liczby falowej wzdłuż tego kierunku. 


Obecność dodatkowego mechanizmu oddziaływającego na przepływ może zmieniać 
i komplikować strukturę turbulencji. Generowane przez parowanie pole sił jest dodatkowym 
źródłem energii kinetycznej turbulencji. Energia ta może być pompowana bezpośrednio 
w małe skale przepływu z ominieciem kaskadowego strumienia energii ([20], [2]). Wpływ 
efektu kondensacji i parowania na strukturę turbulencji jest aktualnie obiektem badań nume- 
tycznych i teoretycznych [21]. Badania te sugerują, że przemiany fazowe mogą wpływać na 
kształt widma energii turbulencji. Brak jednak potwierdzenia eksperymentalnego w jakim 
stopniu te efekty wpływają na globalne charakterystyki struktur przepływu i aglomerację 


kropel w warunkach atmosferycznych. 


Aby głębiej przybliżyć skomplikowaną problematykę omówione zostaną w następnych 
paragrafach zagadnienia takie jak: turbulencja, turbulencja w rzeczywistych chmurach, 


tworzenie się chmur i opadu atmosferycznego, oddziaływanie kropelek chmurowych z prze- 


pływem. 


1.1 Turbulencja 


Fascynacja otaczającym Światem przyrody i prawami 
nim rządzącymi popychała na przestrzeni wieków 
kolejne pokolenia myślicieli do próby opisywania 
i zrozumienia tajemnic natury. Zjawiska zachodzące 
w otaczającej atmosferze i w cieczach intrygowały 
ludzi od tysiącleci. Próby wytłumaczenia i zrozumie- 
nia tych zjawisk można znaleźć u Arystotelesa (Mete- 
orologica). Do naszych czasów przetrwały graficzne 
studia płynącej wody Leonarda da Vinci. Obserwacje 
przepływów turbulentnych były też inspiracją dla 


: r r oe Samp : 
Kartezjusza, który próbował wytłumaczyć sity rządzą- ioe zi) oe COW 
i PRE ; i i Rys. 11: Studium wody  opływającej 
ce ruchem ciał niebieskich przez oddziaływania wirów przeszkodę - Leonardo da Vinci (1508). 


cząstek wypełniających wszechświat. 


Dopiero jednak powstanie dynamiki Newtona i postęp w matematyce pozwoliły na zastoso- 


wanie do badań narzędzi analitycznych. Zastosowanie przybliżenia continuum i rozważania 
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nad wzajemnymi oddziaływaniami infinitezymalnych elementów płynu pozwoliło w końcu na 
zastosowanie dynamiki Newtona do opisu przepływów - równań Naviera - Stokesa (równa- 
nie 1.1). 

Ou 


a, tlt V)a=-V p+vV'i. (1.1) 


Szybko jednak okazało sie, ze te potężne narzędzia da się bezpośrednio wykorzystać do opisu 
tylko pewnej klasy przepływów — przepływów laminarnych. Niestety większość przepływów 
obserwowanych w przyrodzie wykazuje charakter chaotyczny i nie może być analitycznie 


zbadana. 


W drugiej połowie XIX w Osborne Reynolds w czasie eksperymentów z wizualizacją prze- 
pływów w kanale znalazł warunki, w których przepływ z laminarnego przechodził w burzliwy 
ruch turbulentny. Znalazł on bezwymiarową liczbę Re, która określa warunki dla przejścia 


laminarno-turbulentnego [22]: 


Rag. (1.2) 


v 


gdzie Ui L są charakterystycznymi skalami prędkości i długości dla przepływu turbulentne- 
go, natomiast v jest lepkością płynu. Liczba Re jest bezwymiarowym parametrem opisującym 


stosunek pomiędzy bezwładnościowymi i dysypatywnymi właściwościami przepływu. 


Reynolds zaproponował też statystyczną metodę rozwiązania równań ruchu [23]. Zainspi- 
rowany pracami Joule'a potraktował on turbulentne wiry analogicznie jak traktuje się moleku- 
ły w termodynamice statystycznej. Chaotyczna natura przepływów turbulentnych sprawia, że 
z praktycznego punktu widzenia nie interesuje nas ścisłe rozwiązanie przepływu i znalezienie 
chwilowych prędkości, czy trajektorii pojedynczych elementów płynu, podobnie jak nie inte- 
resuje nas rozwiązanie dla położeń pojedynczych molekuł gazu. Istotne są jedynie uśrednione 
wielkości, które charakteryzują przepływ, tak jak w termodynamice istotne są wielkości 
makroskopowe. Reynolds dokonał formalnego podziału przepływu na część średnią i fluktu- 
acje (dekompozycja Reynoldsa). W otrzymanych w ten sposób równaniach ruchu pojawiły się 
dodatkowe naprężenia zwane naprężeniami Reynoldsa. Rozwiązanie tych równań mogło być 
jednak znalezione po dokonaniu pewnych założeń co do tych naprężeń — jest to tzw. problem 


domknięcia układu równań. 
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Rys. 1.2: Turbulentne mieszanie dwóch mas powietrza. Jedna z mas powietrza jest widoczna dzięki 
obecności kropelek wody rozpraszających światło laserowe. Szerokość fotografii odpowiada około 
20cm. Obraz ilustruje różnorodność struktur występujących w przepływie o szerokim zakresie 
rozmiarów. Obraz uzyskany w czasie eksperymentów w komorze chmurowej. 


Z fenomenologicznego punktu widzenia istotną cechą przepływów turbulentnych jest 
występowanie wirowych struktur o różnych wielkościach (skalach przestrzennych, rys. 1.2). 
W pełni rozwinięty przepływ turbulentny zawiera szeroki zakres skal, które są podzielone na 
klasy dużych i małych skal. Duże skale są na ogół związane z wymiarem obszaru czy podsta- 
wowej struktury przepływu. W dużych skalach skumulowana jest większość energii przepły- 
wu. Są one odpowiedzialne za transport pędu, ciepła i masy. Małe skale zawierają przedział 
bezwładnościowy i dysypatywny. Chociaż formalne rozgraniczenie tych skal nie jest trudne, 
to w eksperymentach mamy często problem z niejednoznacznością tego podziału. Źródłem 
energii turbulencji jest energia przepływu głównego, która jest przekazywana największym 
strukturom wirowym występującym w przepływie. Proces ten nazywany jest żargonowo 
„pompowaniem” energii turbulencji. Wskutek oddziaływań inercyjnych największe wiry 
rozpadają się na mniejsze, które z kolei rozpadają się dalej tworząc kolejne generacje coraz 
mniejszych struktur wirowych. Jest to mechanizm przekazywania energii, który prowadzi do 
powstania tak zwanej kaskady energii Richardsona, w której energia kinetyczna jest przekazy- 
wana z większych skal do mniejszych, tylko wskutek bezwładnościowego rozpadu struktur 
wirowych. W rzeczywistych przepływach, wraz ze zmniejszaniem się skali kolejnych genera- 


cji wirów maleje stosunek oddziaływań inercyjnych nad lepkościowymi i inercyjna kaskada 
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energii kończy się na skalach, w których następuje lepka dyssypacja energii kinetycznej. 
Zwyczajowo spektrum skal turbulencji dzieli się na trzy główne części: skalę integralną, 
w której następuje pompowanie energii turbulencji, przedział inercyjny, w którym występuje 
bezwładnościowa kaskada energii i przedział dyssypacyjny, w którym rozpraszana jest energia 
kinetyczna przepływu. W rozwiniętej w pełni, ustalonej turbulencji (stacjonarnej w czasie 
w sensie statystycznym) współczynnik dyssypacji lepkiej <e> jest równy uśrednionemu stru- 
mieniowi energii kinetycznej, która jest pompowana w ruch turbulentny. Jest on również 


równy strumieniowi energii kinetycznej, która jest przekazywana w kaskadzie energii. 


Rozkład energii kinetycznej turbulencji w przestrzeni skal opisuje teoria Kołmogorowa 
[19], która oparta jest na dwóch założeniach. Pierwsze jest takie, że gdy lepkość płynu v jest 
mała to średni współczynnik dyssypacji lepkiej <e> jest niezależny od lepkości. Drugie zało- 
żenie (zwane lokalną izotropią) mówi o tym, że drobnoskalowa turbulencja dla wystarczająco 
dużych liczb Reynoldsa jest statystycznie niezależna od dużych skal i że jest homogeniczna, 
izotropowa i ustalona (tzn. jest statystycznie niezależna od przemieszczeń, obrotów 1 odbić 
układu współrzędnych). Główne wnioski z teorii Kołmogorowa są następujące: statystyczne 
własności najmniejszych skal są określone tylko przez dwie wartości: v i <e>, a dla dużych 
liczb Reynoldsa właściwości statystyczne w skali bezwładnościowej są określone jedynie 
przez <e>. Opierając się na tych założeniach Kołmogorow pokazał na podstawie analizy 
wymiarowej, że funkcja gęstości prawdopodobieństwa unormowanych różnic prędkości 
Au,l(r(e))'*, (gdzie Au,=u(x+r)—u(x) a u jest składową prędkości wzdłuż kierunku x, 
r jest odległością mierzoną wzdłuż x) powinna być uniwersalna, niezależna od przepływu, 
odległości r i liczby Reynoldsa. Można to wyrazić inaczej przez tzw. wzdłużne funkcje struk- 
tury: 

(Au,)=C,(rle))"*, (1.3) 
gdzie n jest liczbą naturalną, a C, - stałą uniwersalną. Podobne równanie można napisać dla 
różnic prędkości z odległością separacji r prostopadłą do kierunku składowej prędkości 


(poprzeczne funkcje struktury). 
Wynika stąd tak zwane prawo ,,4/5” Kołmogorowa: 
(Au))=-Z re), (1.4) 
Odpowiednikiem równania 1.3 dla n=2 w przestrzeni spektralnej jest: 


P(KĘCIE TK. (1.5) 
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gdzie (Kk) jest jednowymiarowym spektrum energii dla liczby falowej x. C jest często 
zwana stałą Kołmogorowa. Na podstawie wyników eksperymentalnych udało się ustalić 


wartość tej stałej (0.5+0.05). 


Zakładając, że właściwości przepływu w małych skalach zależą jedynie od współczynnika 
dyssypacji lepkiej i lepkości płynu, Kołmogorow posługując się analizą wymiarową, znalazł 
charakterystyczne skale dla najmniejszych struktur pola prędkości turbulencji: skalę długości 


(1.6), czasu (1.7) i prędkości (1.8): 


3 1/4 
=| | (1.6) 

v 1/2 
=~], (1.7) 
v=T=|ve]" (1.8) 


Liczba Reynoldsa dla tych skal jest równa jedności: 


1/4 _1/4 3/4 
ovn v Ev 


SĄ p gu Sly (1.9) 


v 


R 


Oznacza to, że dla wirów o rozmiarach porównywalnych ze skalą Kołmogorowa siły inercyj- 


ne są tego samego rzędu co siły lepkościowe. 


Należy zwrócić uwagę, że teoria Kołmogorowa jest fenomenologiczna. Nie wynika ona 
z równań Naviera — Stokesa, ale oparta została na założeniach zainspirowanych obserwacja- 
mi. Powiązanie tej teorii z prawami dynamiki stanowi wciąż wielkie wyzwanie dla współcze- 


snej mechaniki. 
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Rys. 1.3: Kaskada energii zgodnie z teorią Kołmogorowa [22]. Wiry o różnych rozmiarach są 
przedstawione schematycznie jako kontury i ułożone w kierunku malejącego rozmiaru. lo jest skalą 
największych wirów w przepływie, do której pompowana jest energia z ruchu głównego. Skale dla 
kolejnych generacji wirów zmieniają się jak l„=lor", (n=0,1,2,...) gdzie 0<r<1. Liczba wirów na 
jednostkę objętości zmienia się z n jak r”, tak aby kolejne generacje wypełniały całkowicie objętość 
większych wirów. Energia pompowana do największych skal jest następnie przekazywana kolejnym 
generacjom wirów i w końcu dyssypowana w wirach o skali Kotmogorowa n [22]. 


1.2 Turbulencja w chmurach 


Głównym źródłem turbulencji w chmurach konwekcyjnych są przemiany fazowe wody 
w atmosferze. Konwekcja jest wzmagana przez ciepło kondensacji u podstawy chmury. Na 
wierzchołku chmury, wskutek wciągania powietrza suchego do wnętrza chmury, powstają 
obszary o obniżonej wilgotności, w których następuje parowanie kropelek chmurowych. 
Związane z tym lokalne ochłodzenie mas powietrza skutkuje powstaniem przepływu ścinają- 


cego, który jest także źródłem turbulencji [24]. 


Typowa skala największych wirów w przypadku turbulencji w chmurach jest rzędu 10°m. 
Typowe wartości liczby Reynoldsa wynoszą 10° — 10’ przy czym współczynnik dyssypacji 
lepkiej e dla średniej konwekcji w chmurach cumulus jest rzędu 10°m’s” [3]. Turbulencja w 


chmurach w odróżnieniu od większości przepływów występujących w zastosowaniach tech- 
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nicznych charakteryzuje się stosunkowo małą wartością współczynnika e. Z, drugiej strony 


duże skale przepływu skutkują dużymi wartościami liczby Reynoldsa [25]. 


Tabela 1: Charakterystyczne skale Kołmogorowa dla współczyn- 
ników dyssypacji energii typowych dla chmur cumulus[ 26]. 


Współczynnik Skala Czas Prędkość 
dyssypacji energii | Kołmogorowa | Kołmogorowa | Kołmogorowa 

e [m/s] n [m] Tr [S] u, [cm s”] 

10° 4,5x103 1,25 0,35 

107 2,5x10° 0,40 0,63 

10° 1,4x10° 0,13 1,1 

10° 0,8x10° 0,04 2,0 

107 0,45x10° 0,01 3,5 


1.3 Chmury i opad atmosferyczny 


Chmury atmosferyczne to zawiesina wody w postaci ciekłej (kropelki) lub stałej (kryształ- 
ki lodu) unoszonej w powietrzu. Powstają one wskutek procesów kondensacji i resublimacji 
pary wodnej. Do zainicjowania tych procesów konieczne są przesycenia pary wodnej. Do 
takiej sytuacji dochodzi w atmosferze, gdy masa powietrza zostaje ochłodzona poniżej tzw. 
punktu rosy'. Ochłodzenie to jest najczęściej skutkiem adiabatycznego rozprężania wznoszą- 


cej się masy powietrza, mieszania różnych mas powietrza [17], wychładzania radiacyjnego. 


Ciśnienie pary nasyconej zależy nie tylko od temperatury ale też od krzywizny powierzch- 
ni międzyfazowej [24]. Z tego powodu kropelki kondensacyjne nie mogą powstawać w 
atmosferze od zerowego promienia. Ich powstawanie jest zainicjowane przez jądra kondensa- 
cji, czyli drobiny różnego pochodzenia unoszone w powietrzu. Mogą nimi być między innymi 
kryształki soli powstające w czasie falowania morza, pyły powstałe w wyniku erozji czy dzia- 


łalności człowieka, sadza ze spalania biomasy. 


Utworzone na wtrąceniu zarodki kropel w wyniku dyfuzji pary wodnej i jej kondensacji 
powiększają swój rozmiar. Proces ten staje się coraz mniej wydajny wraz z powiększaniem 
się kropelki co wynika ze stosunku promienia kropelki do jej objętości (~r’). Dla wzrostu 


kropelek powyżej 36um znaczenia nabierają zderzenia między kropelkami i koalescencja 


I Temperatura punktu rosy - temperatura, w której przy danym ciśnieniu i wilgotności para wodna zawarta w 
powietrzu staje się nasycona, a poniżej tej temperatury staje się przesycona i skrapla się. 
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kropel ([24], [16], [27]). Ten mechanizm jest głównie odpowiedzialny za powstawanie opadu 


atmosferycznego. 


Do powstania jednej kropli opadowej potrzebna jest jednak bardzo długa sekwencja 
zderzeń. Istotną rolę w tym procesie odgrywa bezwładność kropelek i ich oddziaływania 


z przepływem, a szczególnie z turbulentnym polem prędkości [14]. 


Problem ten nie jest nadal rozwiązany i zrozumienie jak rozkład wielkości kropelek chmu- 
rowych zmienia się w czasie i w przestrzeni jest jednym z największych wyzwań we współ- 


czesnej fizyce chmur. 


1.4 Oddziaływania kropelek z przepływem 


Kropelki chmurowe nie są pasywne, oddziaływają z otaczającym powietrzem na kilka 


sposobów: 
e oddziaływania hydrodynamiczne, 
e oddziaływania termodynamiczne, 
e wymiana masy. 


Kropelki chmurowe, które zawierają dodatkowe substancje chemiczne mogą oddziaływać 
również z otoczeniem poprzez reakcje chemiczne. Problem ten nie jest jednak rozważany 


w tej pracy. 


Oddziaływania hydrodynamiczne polegają na wymianie energii mechanicznej pomiędzy 
płynem a kropelką. W dynamicznym przepływie następuje ciągła wymiana energii kinetycz- 
nej i potencjalnej pomiędzy powietrzem i unoszoną w nim kropelką. Różnica gęstości pomię- 
dzy dwiema fazami skutkuje różną bezwładnością kropelki i otaczającego ją medium, która 
powoduje, że kropelki poruszają się po trajektoriach nieco innych niż otaczający je płyn. 
W ten sposób mogą one pobierać energię kinetyczną przepływu w miejscach, w których są 


przyspieszane i oddawać ją w innym miejscu, w którym wyhamowują. 


W przypadku rozważanych tu kropelek o rozmiarach rzędu 10 mikrometrów, różnice 
między ruchem kropelki i otaczającego ją płynu jest niewielka. Jest to spowodowane relatyw- 
nie małym stosunkiem sił objętościowych do wypadkowej siły działającej na powierzchnię 
kropelki. Podczas gdy siły objętościowe (bezwładność i grawitacja) dążą do zmiany prędkości 
kropelki względem płynu, siły oporu działające na powierzchnię wyrównują prędkość kropel- 
ki do prędkości otaczającego ją płynu. Siły objętościowe zależą od promienia kropli jak r” 
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natomiast wypadkowa siła powierzchniowa (siła Stokesa) jak r. Zatem stosunek sił objęto- 
ściowych do sił powierzchniowych zmienia się wraz z wielkością jak 7°. Zależność ta została 


ilościowo wyrażona przez tzw. czas Stokesa z, [26] (czas relaksacji dla sedymentacji): 


2 
m 
tT, =—— ETEN (1.10) 
P? 6nrRu u 


gdzie m, jest masą kropelki, R jej promieniem, u lepkością dynamiczną płynu, pw gęstością 
wody. 


Jest to czas, po którym opadająca swobodnie w polu grawitacyjnym kropelka osiągnie 


około 63% prędkości granicznej. 


Tabela 2: Skale charakterystyczne dla wybranych wielkości kro- 
pelek wody w powietrzu [26]. 


A Liczba Prędkość Czas 
R [m] Reynoldsa końcowa Stokesa 
Rep V; [cm/s] Tp [S] 

5x106 0,9x10* 0,3 3x10" 
10x10° 0,7x10° 1,2 11x104 
15x10% 0,2x10"! 2] 26x104 
20x10% 0,6x10" 4,8 46 x10-4 


Małe prędkości kropelek względem przepływu skutkują małymi liczbami Reynoldsa (tabe- 
la 2). Można przyjąć, że dla kropelek wody w powietrzu obowiązuje prawo Stokesa dla siły 
oporu: 


F=6ru R(u-v), (1.11) 


gdzie u jest lepkością dynamiczną powietrza, R promieniem kropelki, u prędkością prze- 


pływu a v prędkością kropelki. 


Drugie prawo dynamiki Newtona dla sfery poruszającej się z prędkością v w jednorodnym, 


ale zmiennym w czasie przepływie o prędkości u wygląda następująco [28], [3]: 


dv 


u(t')-v(t') 
= ; 
Pa d dt t—t' 


I > ( 
=6ruR(u-v)+>p,V,(ui—V)+6R ym — = di 
uR( ) 274 al ) pe] a (1.12) 
+p,V .g+p,V,(ii—g). 
gdzie u jest lepkością płynu, py; jego gęstością, pa to gęstość cząsteczki, a V„=4/3nR* jest obje- 


tością cząsteczki. Po prawej stronie równania 1.12 występują następujące człony: siła Stoke- 
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sa, dodatkowa siła masowa spowodowana przyspieszeniami w otaczającym płynie, siła 
Basseta spowodowana dyfuzją wirowości powstającej podczas przyspieszania cząsteczki, siła 
grawitacji. Ostatni człon to poprawka na ścinanie pola prędkości wskutek oddziaływania na 


cząsteczkę i siłę wyporu. 


Równanie to powstało wskutek uproszczeń takich jak między innymi: zaniedbanie zakrzy- 
wienia pola prędkości, zaniedbanie poślizgu, oddziaływań powierzchniowych i oddziaływań 


z innymi kropelkami 


Oddziaływania hydrodynamiczne kropelek z otaczającym je przepływem odbywa się 
w dwóch kierunkach: przepływ powietrza powoduje zmianę trajektorii kropelek, kropelki 
powodują lokalne zmiany pola prędkości. Ze względu na małe liczby Reynoldsa pojedyncze 
kropelki nie produkują turbulentnych wirów. Rozpatrywane przepływy w chmurach są 
w ogólności turbulentne, jednak skala Kołmogorowa jest dwa, trzy rzędy wielkości większa 
od rozmiarów kropelek (porównaj tabele 1 i 2). To sprawia, że przepływ w otoczeniu kropelki 
jest wystarczająco gładki, aby kropelka nie „zauważała” przestrzennych zmian pola prędko- 
ści. W skalach oddziaływań pojedynczej kropelki z przepływem (poniżej skali Kotmogorowa) 
przepływ nie jest turbulentny. Stosowanie prawa Stokesa i równania 1.12 do opisu ruchu 


kropelek chmurowych jest więc w tym przypadku uzasadnione. 


Ponieważ wielkość kropelek jest mała w porównaniu z typowymi wartościami skali 
Kołmogorowa w chmurach, oznacza to, że kropelki nie oddziaływają bezpośrednio z wirami 
turbulentnymi, bo są mniejsze co najmniej dwa rzędy wielkości od najmniejszych wirów 
występujących w przepływie. Wynika z tego również, że wiry turbulentne i związanie z nimi 


gradienty prędkości nie powodują deformacji powierzchni kropelek i ich rozrywania [1]. 


Oddziaływania hydrodynamiczne pomiędzy przepływem a kropelką mogą być ilościowo 


opisane przez tzw. liczbę Stokesa: 
S=—. (1.13) 


Liczba ta wyraża stosunek czasu Stokesa t, do charakterystycznej skali czasu zmienności 
przepływu Tr co daje liczbowe oszacowanie szybkości z jaką cząsteczka dostosowuje się do 
zmian pola przepływu. Mała liczba Stokesa odpowiada sytuacji, gdy kropelka porusza się 
niemalże jednocześnie z płynem. Duża wartość tej liczby odpowiada sytuacji, gdy pole pręd- 
kości płynu w nikły sposób wpływa na prędkość cząsteczek. Spodziewamy się największych 
oddziaływań pomiędzy przepływem i kropelką dla liczby Stokesa równej jedności. 


1.4 Oddziaływania kropelek z przepływem 


W przypadku przepływu turbulentnego charakterystyczną skalą czasu dla przepływu może 
być skala Kołmogorowa t», ponieważ w tej skali spodziewamy się maksymalnej wirowości 
czyli maksymalnych zmian pola prędkości. Tabela 3 zawiera liczby Stokesa obliczone dla 


wybranych promieni kropel i czasów Stokesa (porównaj tabela 2). 


Tabela 3: Liczba Stokesa dla wybranych promieni kropelek 


i wybranych skal czasu Kołmogorowa. 


R [m] 
5x10*m 


Tr [S] 


10x10 m 


15x10% m 


20x10% m 


1,25 s 0,0002 | 0,0009 0,0021 0,0036 
0,40 s 0,0007 | 0,0028 0,0064 0,0114 
0,13 s 0,0022 | 0,0088 0,0197 0,0351 
0,04 s 0,0071 0,0285 0,0641 0,1140 
0,01 s 0,0285 0,1140 0,2564 0,4558 


Jak można zauważyć dla typowych skal turbulencji występujących w chmurach (por. Tabe- 
la 1) liczba Stokesa jest mniejsza od jedności. Możemy zatem założyć, że tory kropelek 
chmurowych wiernie odtwarzają lokalne struktury pola prędkości. Nie należy jednak zapomi- 


nać o wpływie wymienionych wcześniej sił wypornościowych, modyfikujących ruch struktur 


wirowych. 


2 Diagram mieszania — analiza termodynamicznych 
efektów procesu mieszania 


Jak to zostało opisane w poprzednim rozdziale, obecność wody w atmosferze czy to 
w stanie ciekłym czy gazowym prowadzi do powstania mechanizmów, które mogą zmieniać 
przepływ. Oprócz oddziaływań hydrodynamicznych pojawiają się również efekty związane 
z parowaniem, kondensacją i dyfuzją pary wodnej. Prowadzi to między innymi do zmian 
gęstości, co jest bardzo istotne w przepływie turbulentnym, ponieważ fluktuacje gęstości 


prowadzą do wygenerowania fluktuacji sił wyporu, które oddziaływają na przepływ. 


W obecności turbulencji, masy powietrza charakteryzujące się różnymi własnościami 
mieszają się w różnych proporcjach, tworząc skomplikowane przestrzenne struktury (patrz 
rys.1.2). W niniejszym rozdziale opisany zostanie prosty model mieszania dwóch różnych 
mas powietrza w ustalonym ciśnieniu. Do opisu zmian gęstości użyta zostanie temperatura 


gęstościowa T, zdefiniowana jako[29]: 


R, 
TORTNE = Aw 


a 


l (2.1) 


gdzie 7 jest fizyczną temperaturą powietrza, R,, R, stałymi gazowymi pary wodnej i suchego 
powietrza (powietrza o wilgotności względnej 0%), X” stosunkiem zmieszania pary wodnej, 


w™ masowa zawartość wody kropelkowej. 


Tak zdefiniowana temperatura gęstościowa jest proporcjonalna do odwrotności gęstości: 


1 
PR (2.2) 


p 
Uwzględnia ona dwie poprawki gęstości paczki powietrza atmosferycznego, spowodowane 


obecnością wody. Pierwsza poprawka — 1+(1-7)X uwzględnia wpływ pary wodnej na 


gęstość powietrza, druga - w - uwzględnia wpływ wody kropelkowej zawartej w analizowanej 


objętości. 


I R= R,/ u,,R=R al U, , gdzie R g Jest uniwersalną stałą gazową, a H,,, Ha masami molowymi pary 
wodnej i suchego powietrza. Pozwala to na zapisanie równania stanu składnika n ( p, V =N „RT ) przy 
użyciu parametrów intensywnych ( p,=p,R,„T'). 

II Dowolną objętość powietrza można traktować jako mieszaninę powietrza suchego i pary wodnej. Stosunek 


zmieszania X jest ilorazem masy pary wodnej m, do masy suchego powietrza mą: X=m,/mg. 
HI Jeżeli m jest masą danej objętości powietrza, a m,, masą wody kropelkowej w tej objętości to w=m,,/m. 
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Rozważmy dwie różne masy powietrza o ustalonych parametrach termodynamicznych. 


Nazwijmy je odpowiednio powietrzem czystym i powietrzem chmurowym. 


Powietrze czyste niech będzie powietrzem pozbawionym kropelek wody. Może zawierać 
parę wodną, stanowi zatem mieszaninę suchego powietrza i pary wodnej. Skład procentowy 
obu składników mieszaniny jest jednoznacznie określony przez stosunek zmieszania (por 


rozdz. 9.1). 


Powietrze chmurowe jest mieszaniną suchego powietrza i pary wodnej, dodatkowo unoszą 
się w nim kropelki wody. Masa ciekłej wody w postaci kropelek znajdujących się w danej 
objętość powietrza chmurowego jest określona przez masową zawartość wody kropelkowej 
w. Kropelki są w równowadze termodynamicznej z powietrzem, co oprócz warunku równości 
temperatury narzuca warunek na stosunek zmieszania pary wodnej. Zawartość pary wodnej 
musi odpowiadać stanowi nasycenia dla danej temperatury i ciśnienia. To zapewnia, że 


kropelki nie parują w takich warunkach. 


Zastanówmy się jakie parametry termodynamiczne ustalą się po jednorodnym zmieszaniu 
tych obu porcji powietrza dla ustalonego stosunku zmieszania powietrza chmurowego K, 
który definiuje się jako stosunek masy powietrza chmurowego do sumy obu porcji powietrza 
(chmurowego i czystego). Wartość k=0 odpowiada przypadkowi, gdy mamy tylko czyste 


powietrze, wartość k=1 odpowiada przypadkowi, gdy mamy tylko powietrze chmurowe. 


Zbadajmy jaka temperatura gęstościowa ustali się w rozważanym przez nas procesie 
mieszania. Sporządzony diagram mieszania dla parametrów eksperymentalnych (T=298K) i 
różnych wilgotności względnych czystego powietrza przedstawiony jest na rys.2.1. Przedsta- 
wia on zależność pomiędzy temperaturą gęstościową a proporcjami zmieszania powietrza 
chmurowego 7,(k). Możemy zauważyć, że dla stałej temperatury powietrza czystego tempera- 
tura gęstościowa 7, dla różnych wilgotności zmienia się. Jest to spowodowane tym, że zmia- 
na wilgotności powoduje zmianę zawartości pary wodnej (por. 2.1). W przeprowadzonych 
eksperymentach, z uwagi na stałe warunki termodynamiczne temperatura gestościowa 7, była 


stała. 
Wyraźnie widać, że wszystkie wykresy na diagramie mieszania składają się z dwóch części 
i charakteryzują się skokiem pochodnej w punkcie łączącym obie części. Pierwsza część 


wykresu (dla małych wartości k) jest we wszystkich przypadkach malejąca, natomiast nachy- 


I Jeżeli masę porcji powietrza chmurowego oznaczyć przez m,, a przez m» masę porcji powietrza czystego, to 


m, 


proporcje zmieszania będą siĘ Wyrazac wzorem: k = m,+m, 
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lenie drugiej części wykresu (dla dużych wartości k) zmienia się wraz z wilgotnością. Obie 
części wykresu odpowiadają dwóm różnym reżimom zależnym od początkowych parametrów 


termodynamicznych. 


T, [K] 


0 0.2 0.4 0.6 0.8 1 


Rys. 2.1: Diagram mieszania — zależność temperatury gęstościowej Tp od proporcji zmieszania k 
chmurowego powietrza z powietrzem czystym. Krzywe odpowiadają różnym wartościom wilgotności 
względnej czystego powietrza przed zmieszaniem (od 20% do 80%). 


Pierwszy reżim odpowiada sytuacji, gdy mieszana jest niewielka porcja powietrza chmuro- 
wego, więc jest na tyle mało kropelek, że wszystkie parują i równowaga termodynamiczna 


ustala się dla wilgotności mniejszej od RH=100%. 


Drugi reżim występuje, gdy wody kropelkowej jest na tyle dużo, że stan równowagi termo- 
dynamicznej jest osiągany zanim odparują wszystkie kropelki, wtedy w rezultacie zostaje 
osiągnięty stan nasycenia (wilgotność względna równa 100%). Część kropelek, które nie 
odparowały zupełnie pozostaje nadal w formie zawiesiny kropel będących w równowadze 


termodynamicznej z otaczającym powietrzem. 


Punkt, w którym następuje skok pochodnej odpowiada sytuacji pośredniej, gdy po odparo- 


waniu całej wody kropelkowej zostaje osiągnięte nasycenie. 


Omawiany diagram ilustruje jak skomplikowany jest proces rozważanego mieszania 
powietrza chmurowego. Można sobie wyobrazić, że w wyniku tego procesu masy powietrza 
o różnej zawartości wody mieszają się niejednorodnie (patrz rys. 1.2). Powoduje to występo- 
wanie obok siebie porcji powietrza różniących się historią mieszania, co skutkuje powstaniem 
różnicy wyporu pomiędzy sąsiednimi obszarami i pojawieniem się fluktuacji siły wyporu w 


mieszającej się objętości. 


Zarówno dla chmury atmosferycznej jak i w warunkach laboratoryjnych nie znamy rozkła- 


du gęstości prawdopodobieństwa dla stosunku zmieszania k. Zakładając najprostszy przypa- 
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dek, to znaczy rozkład płaski, w którym każda wartość k jest równie prawdopodobna i wystę- 
puje równie często, możemy przeanalizować statystykę temperatury gęstościowej Tp. Jest to 
zapewne bardzo duże uproszczenie, ale przy braku danych przynajmniej pozwoli nam na 
jakościową dyskusję diagramu mieszania. Uśredniając wykres Tọ można oszacować jak zmie- 
nia się średnią gęstość mieszającej się chmury. Rys. 2.2 przedstawia różnicę pomiędzy Tp dla 
powietrza czystego a średnią wartością 7, dla całego diagramu: 
AT,=T,(0)=(T,(k)). (2.3) 

Operacja <> oznacza tutaj uśrednianie po wszystkich wartościach proporcji zmieszania 
ke|0,1|. Możemy przypuszczać, że wykres ten pokazuje jakościowo, jak zmienia się średnia 
różnica wyporności pomiędzy powietrzem w opadającej chmurowej strudze a otaczającym 
czystym powietrzem w zależności od wilgotności względnej czystego powietrza. Zależność ta 
sugeruje, że dla mniejszych wilgotności powinna pojawiać się średnio większa różnica wypo- 


ru co powinno skutkować większą średnią prędkością opadania powietrza chmurowego. 


i i i i i i 
10 20 30 40 50 60 70 80 90 100 
wilgotnosc wzgledna [%] 


Rys. 2.2: Średnia różnica temperatury gęstościowej między powietrzem suchym a zmieszanym 
w zależności od wilgotności. 


Lokalne fluktuacje siły wyporu opisują fluktuacje temperatury gęstościowej. Przyjmiemy, 
że wyodrębnione porcje powietrza bezpośrednio ze sobą sąsiadujące będą miały podobną 
historię mieszania (inaczej mówiąc powstały przez wymieszanie dwóch różnych mas powie- 


trza z podobnymi proporcjami zmieszania k). Lokalne zmiany siły wyporu zależą od pochod- 


. pzd . . : . dT, : 3 . 
nej temperatury gęstościowej Tp, po proporcjach zmieszania k - =~. Z diagramu mieszania 


wynika, że największe pochodne występują dla pierwszej części wykresu i zakres zmian k jest 
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większy dla mniejszych wilgotności. Ilustruje to moduł pochodnej dla diagramu mieszania 


pokazany na rys.2.3. 


RH=30% 
RH=40% 
RH=50% 
= = = RH=60% 
RH=70%|- 
RH=80% 


Modul pochodnej diagramu mieszania [K] 


m R MZ m m in m m in mm m im mm I 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
k 


Rys. 2.3: Moduł pochodnej diagramów mieszania w funkcji stosunku zmieszania k. dla różnych 
wilgotności względnych RH. 
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Rys. 2.4: Diagram mieszania dla wilgotności wzglednej w okolicy RH=66%. 

Wynika z tego, że największe różnice wyporu powstają gdy mieszane jest powietrze czyste 
z małą proporcją powietrza chmurowego, można się więc spodziewać silnych efektów na 
brzegu chmury. Efekty te będą silniejsze, gdy różnica wilgotności pomiędzy obiema masami 
będzie duża, wtedy wydłuża się obszar dużych pochodnych Tp. Różnice modułu pochodnej 
w zależności od wilgotności są mało istotne dla pierwszej części wykresu, natomiast dla 
drugiej części wykresu pochodna zmienia się w sposób niemonotoniczny. Na diagramie 


mieszania (rys. 2.1) można zauważyć, że wraz ze zmianą wilgotności zmienia się znak 
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pochodnej dla drugiej części wykresu. Najmniejsze wartości modułu pochodnej występują dla 


wilgotności około 66% (por. rys. 2.4). 


0.4 


RH=30% 
RH=40% 
RH=50%| | 
= = = RH=60% 
RH=70%] 
RH=80% 


czestosc wystepowania 


Rys. 2.5: Rozkłady temperatury gęstościowej dla różnych wilgotności względnych RH sporządzone na 
podstawie diagramu mieszania. 


W analizie fluktuacji wyporu trzeba wziąć też pod uwagę stykające się cząstki powietrza, 
które mają inną historię mieszania. W tym przypadku należy przeanalizować histogram 75. 
Na pierwszy rzut oka można zauważyć na rys. 2.1, że krzywe dla wilgotności względnej 
RH=20% i RH=80% różnią się istotnie. Wyraźnie widać na przykład, że zakres zmian warto- 
ści Tp dla mniejszej wilgotności jest o wiele szerszy niż zakres zmian wartości dla większej 
wilgotności. Tendencję tę można zauważyć również na histogramach wyznaczonych dla 
rozkładu Tp. (rys. 2.5). Przedstawia on unormowane do jedności rozkłady 7, dla różnych 
wilgotności. Rozkłady różnią się w znaczny sposób między sobą zarówno szerokością jak 
i kształtem. Dla wilgotności w okolicach RH=66% rozkład staje się najwęższy. Szerokość 
rozkładów została oszacowana przez odchylenie standardowe (rys. 2.6). Wynika z tego, że 


szerokość rozkładu maleje wraz z wilgotnością do RH~66% następnie zaczyna rosnąć. 


Podsumowując, w przedstawionych powyżej rozważaniach zbadano na podstawie diagra- 
mów mieszania proces mieszania powietrza chmurowego z powietrzem suchym. Pomimo 
grubych uproszczeń (założony rozkład płaski proporcji zmieszania k) można wyciągnąć 
pewne jakościowe wnioski. Wynika z nich, że prędkość opadania paczki mieszającego się 
powietrza powinna rosnąć wraz ze wzrostem różnicy wilgotności pomiędzy dwiema masami 
powietrza. Natomiast szerokość rozkładu fluktuacji gęstości zależy od wilgotności w sposób 
niemonotoniczny i może się zmieniać w szerokim zakresie w zależności od lokalnych warun- 


ków termodynamicznych. 
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odchylenie standardowe dla rozkladu Tp [K] 


0 10 20 30 40 50 60 70 80 90 100 
wilgotnosc wzgledna [%] 


Rys. 2.6: Wykres zależności odchylenia standardowego dla rozkładu temperatury gęstościowej 
w zależności od wilgotności względnej. 
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3 Anemometria obrazowa - PIV 


3.1 Opis metody 


Podstawową techniką badania pól prędkości zastosowaną w tej pracy jest anemometria 
obrazowa - PIV (z angielskiego — Particle Image Velocimetry) [30]. Technika ta, dzięki opar- 
ciu na analizie obrazów, posiada niewątpliwy atut jakim jest nieinwazyjność, przy czym 
dostarcza bogatych danych w postaci pól prędkości czy nawet pól przyspieszeń. Nic dziwne- 
go, że znalazła zastosowanie w różnorodnych dziedzinach, od tych bardziej oczywistych 


zastosowań w technice przy projektowaniu elementów pojazdów po biologię i medycynę. 


Różnorodność zastosowań i postęp technologiczny sprzyjały wykształceniu się wielu tech- 


nik pochodnych jak na przykład 3D-PIV, Stereo-PIV czy Micro-PIV. 


Eksperymenty z wykorzystaniem meto- 
dy PIV polegają na utrwalaniu chwilo- 
wych obrazów wizualizowanego przepły- 
wu. Aby ruch płynu stał się widoczny 
rozprowadza się w nim cząsteczki znacz- 
nika na tyle małe, aby nie zmieniały prze- 
pływu i jednocześnie poruszały się razem 


z nim. Badany przepływ oświetlany jest 


silnym źródłem światła, które rozpraszane Rys. 3.1: Uproszczony schemat układu eksperymentalnego 


jest przez cząsteczki. Światło jest formo- P/V! - źródło światła, 2 — układ formujący płaszczyznę 
świetlną 3 — płaszczyzna świetlna podświetlająca 
wane przy użyciu soczewki cylindrycznej fragment przepływu. 
w płaską, cienką wiązkę (nóż świetlny). Dzięki temu otrzymywany jest prawie dwuwymiaro- 
wy przekrój przez badany przepływ, który jest rejestrowany przez kamerę ustawioną prosto- 
padle do płaszczyzny przekroju. W czasie eksperymentów PIV pozyskiwane są pary fotogra- 
fii, przy czym dwa zdjęcia w parze otrzymywane są w różnych chwilach czasu. Współcześnie 
jako źródła światła używa się najczęściej laserów impulsowych. Oświetlenie przepływu krót- 
kimi, ale intensywnymi impulsami światła pozwala na rejestrację nierozmytych obrazów 


cząstek znacznika przy dużych prędkościach przepływu. Zastosowanie kamer cyfrowych 


pozwala na automatyczną akwizycję dużych liczb fotografii, które są gotowe do dalszej anali- 
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zy przy pomocy komputerowych programów PIV. Fotografie takie pozwalają na utrwalenie 


chwilowego stanu położeń cząsteczek znacznika. 


Układ obrazujący 
kamery 


Smuga światła o 
grubości AZ 
X oświetlająca fragment 


= przepływu 
y 


Płaszczyzna fotografii 
Rys. 3.2: Schemat rzutowania na płaszczyznę fotografii. Trójwymiarowy element objętości jest 
rzutowany przez układ obrazujący kamery w płaski fragment fotografii. 

Przeanalizujemy metodykę wyznaczania pola prędkości na podstawie obrazów cząstek. 
Niech podświetlony fragment przepływu zawiera zbiór N cząsteczek. Chwilowy rozkład prze- 
strzenny cząsteczek jest opisany przez zbiór I: 

R, 
R, 
PEP. (3.1) 


R 


N 
gdzie wektor R; = [X,, Y,, Zi] jest wektorem położenia i-tej cząsteczki w trójwymiarowej prze- 
strzeni przepływu. Jeśli zaniedba się grubość wiązki świetlnej to: 

V Z,=const | (3.2) 
więc położenie widocznych cząsteczek opisywane jest tylko przez dwie składowe Xi, i Y;. 


Rozproszone na cząsteczkach znacznika Światło rzutowane jest przez układ optyczny 
kamery na matrycę CCD (rys. 3.2), która rejestruje obraz i wysyła go do urządzenia przecho- 


wującego dane w postaci dwuwymiarowej cyfrowej fotografii. Cząsteczkom znacznika 
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w przestrzeni przepływu odpowiadają obrazy tych cząsteczek na dwuwymiarowej płaszczyź- 
nie fotografii. Zbiorowi / odpowiada zbiór y: 

ry 
F3 


y=("|, (3.3) 


gdzie r; = [x, yj] jest położeniem j-tego obrazu cząsteczki na płaszczyźnie fotografii 
a n liczebnością zbioru. Każdej cząsteczce przyporządkowany jest obraz cząsteczki na płasz- 
czyźnie fotografii. W ogólności nie jest to jednak odwzorowanie różnowartościowe i jedne- 
mu obrazowi może odpowiadać kilka cząsteczek w rzeczywistej przestrzeni. Jest to spowodo- 
wane tym, że kilka cząsteczek w trójwymiarowej przestrzeni może być rzutowanych w to 
samo miejsce na płaszczyźnie fotografii (cząsteczki mogą się zasłaniać), więc wielkości zbio- 
rów spełniają następującą relację: 
n<N. (3.4) 
W eksperymentach z użyciem techniki PIV staramy się używać na tyle cienkich smug 


świetlnych aby efekty przestrzenne były pomijalne. Dąży się więc do tego aby: 
N-n 


«x1. (3.5) 


Jeśli spełnione są te warunki to oznacza to, że prawie każda cząstka znacznika ma swój 


odpowiednik na płaszczyźnie fotografii. 


Współrzędne na płaszczyźnie fotografii związane są ze współrzędnymi rzeczywistymi 


przekształceniem liniowym: 


X=—,Y=.. (3.6) 


sal 
M? M 

Współczynnik M jest współczynnikiem skalowania, który jest znajdowany empirycznie 
w procesie kalibracji układu pomiarowego. 

Pojedyncza fotografia zawiera więc informacje o chwilowym stanie układu, czyli chwilo- 


wym rozkładzie przestrzennym cząsteczek znacznika. 


I W rzeczywistości mamy do czynienia z rozkładem jasności obrazu cząstki a jej położenie wyznacza punkt 
centralny tego rozkładu 
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W eksperymentach PIV rejestrowane są pary fotografii'. Pierwsza fotografia wykonywana 
jest w chwili to, następna w t=t)+At. Podobnie jak to było zrobione wcześniej, dla każdej 
z chwil czasowych można zdefiniować zbiory położeń cząsteczek znacznika dla każdego 
z obrazów: 7 (t9) i T(t). W praktyce, jeżeli At40 to oba zbiory nie dotyczą tego samego zbioru 
cząsteczek. W przedziale czasowym At część cząsteczek ucieknie poza kadr kamery, a jedno- 
cześnie pojawią się inne, które nie były widoczne wcześniej. To samo dotyczy również płasz- 
czyzny świetlnej, w którą wpadają i wypadają cząsteczki. W zależności od długości At, grubo- 
ści płaszczyzny świetlnej 1 od charakteru przepływu będzie zmieniała się liczba cząsteczek, 
które widoczne są dla obu chwil czasowych. Zbiory 7 (ty) i I (t,) zapiszemy w postaci sumy 
podzbiorów: 


T(t,) 
T(t,) 


Zbiory oznaczone indeksem w zawierają położenia tego samego zbioru cząsteczek 


=T,(eJUT, (1). ee 
w dwóch różnych chwilach czasu. Maja one zatem tę samą liczbę elementów N,. Zbiory 
oznaczone indeksem s zawierają położenia cząsteczek, które widoczne są tylko w jednej 


chwili czasowej. 


Podobnie jak to zostało pokazane wcześniej, zbiorom położeń cząsteczek w rzeczywistej 
przestrzeni odpowiadają zbiory położeń obrazów tych cząsteczek na płaszczyźnie fotografii: 


Y (4%)=Y,(%0UY,(t0). 


y(t,)=y,(t,)Uy,(t). (3.8) 


Para fotografii zawiera więc informację o położeniach tych samych cząsteczek znacznika 
w różnych chwilach czasu — zbiory z indeksem w i informację o położeniu przypadkowych 


cząsteczek — zbiory z indeksem s. 


Jeżeli wielkość zbiorów y, jest mniejsza od wielkości zbiorów y, to dwie fotografie zawie- 
rają informację o przemieszczeniach cząsteczek, którą można odzyskać przez przetwarzanie 
tych fotografii technikami PIV. Praktycznie wymagane jest aby przynajmniej 60% cząstek 


znajdowało się na obu obrazach. 


I Istnieją również metody PIV wykorzystujące wiele obrazów, ale w tej pracy nie były stosowane. 
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3.2 Opis programu PIV-Kor zastosowanego w badaniach 


Wstępne badania przepływu w modelu chmury wskazały, że konieczne jest stworzenie od 
podstaw własnego oprogramowania PIV. Eksperymenty w komorze chmurowej różnią się 
istotnie od klasycznych eksperymentów z wykorzystaniem PIV, gdzie do przepływu wprowa- 
dza się specjalnie dobrane cząstki znacznika w ten sposób, żeby równomiernie wypełniały 
objętość pomiarową. W przypadku procesu mieszania w komorze chmurowej kropelki, które 
pełnia rolę znaczników, nie wypełniają przestrzeni równomiernie a silne gradienty ich 
koncentracji mogą stanowić poważną trudność w stosowaniu techniki PIV. Dostępne oprogra- 
mowanie PIV zastosowane do fotografii z komory chmurowej nie dawało zadowalających 
i wiarygodnych wyników. Autor zapoczątkował więc pracę nad programem PIV-Kor, który 
nadal jest rozwijany. Oprogramowanie mimo braku przyjaznego interfejsu okazało się bardzo 
użyteczne w przygotowaniu pracy, napisane w środowisku MATLAB umożliwia dosyć łatwe 
wprowadzanie modyfikacji, które pozwalają na wprowadzanie w nim ulepszeń i dostosowy- 
wanie do konkretnych zastosowań. Dzięki temu program PIV-Kor mógł być stosowany 
również w innych badaniach (patrz np. [1]). Samodzielna praca nad realizacją algorytmu PIV 
pozwoliła również autorowi na dogłębną analizę stosowanej w badaniach techniki pomiaro- 
wej i bardziej krytyczne do niej podejście. W ostatnim czasie dostępne stało się również 
nowoczesne oprogramowanie firmy ILA GmbH — VidPIV. Porównanie działania obu progra- 
mów (PIV-Kor i VidPIV) pokazało, że dają one niemal identyczne wyniki (patrz rozdz. 9.3), 
co potwierdza poprawność wybranego podejścia. Z, drugiej strony uwiarygodniło to rezultaty 
oprogramowania VidPIV, stosowanego w tej pracy do pomiarów wykonanych techniką 


Stereo-PIV (patrz rozdz. 5.1). 


W tym rozdziale skupimy się na omówieniu zasady działania programu PIV-Kor. Bardziej 


szczegółowe omówienie algorytmu programu znajduje się w dodatkach (rozdz. 9.2). 


Głównym elementem programów stosujących technikę PIV jest algorytm wyszukiwania 
lokalnych przemieszczeń, to znaczy przemieszczeń pojedynczych cząsteczek znacznika lub 
pojedynczych pikseli obrazu. W większości algorytmów PIV stosuje się raczej to drugie 
podejście i szuka się przemieszczeń małych elementów obrazów. W ten sposób punktowi 
obrazu można przypisać przemieszczenie wyznaczone dla fragmentu będącego otoczeniem 
tego punktu. Pomiar prędkości na podstawie dwóch fotografii przepływu jest możliwy, gdy 
przemieszczenia cząsteczek są niewielkie. Inaczej sprawiałoby to poważne błędy interpreta- 


cyjne. Dzięki temu można założyć, że grupa cząsteczek znajdujących się na fotografiach nie 
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ulega znacznym deformacjom, jest tylko przesunięta. To założenie pozwala na sprowadzenie 


problemu znajdowania przemieszczeń do problemu porównywania obrazów. 


Idea programów PIV bazuje na algorytmach przetwarzających obrazy. Podstawową opera- 
cją takich algorytmów jest wyszukiwanie podobieństw dla kolejnych obrazów przepływu. Jest 
to konieczne, gdy chcemy znaleźć położenia tych samych obiektów na obu fotografiach. 
Obiektowi w rzeczywistej przestrzeni odpowiada obraz tego obiektu na płaszczyźnie fotogra- 
fii. Dwie fotografie na których zarejestrowano ruch obiektu zawierają obrazy tego obiektu 
uchwycone w różnych położeniach. Obrazy obiektu zarejestrowane w różnych chwilach czasu 
są podobne ale praktycznie nigdy nie są identyczne (jest to związane ze specyfiką samego 
zjawiska, w którym sam obiekt może podlegać deformacjom, obrotom, może to wynikać 
również z przyczyn technicznych takich jak różnica w oświetleniu, szumy układu obrazujące- 
go itp.). Do porównywania obrazów stosuje się różnorodne metody opisane w literaturze [30]. 
Program stworzony na potrzeby niniejszej pracy oparty jest na współczynniku kowariancji 


lub ściślej na współczynniku korelacji liniowej. 


Cyfrowe obrazy można traktować jak tabele, czy macierze. Indeksy poszczególnych 
elementów macierzy są związane liniowym przekształceniem ze współrzędnymi na płasz- 
czyźnie fotografii i będą tutaj utożsamiane jako współrzędne - odpowiednio indeks i — współ- 
rzędna pozioma, j - pionowa. Wartości poszczególnych elementów macierzy odpowiadają 
poziomowi jasności w danym punkcie fotografii. Jasność przybiera całkowite wartości dodat- 
nie z zakresu określonego przez właściwości kamery i liczbę bitów przetwornika analogowo- 


-cyfrowego. 


Niech $7, S2 oznaczają odpowiednio macierze odpowiadające cyfrowym obrazom o tych 
samych rozmiarach MxN. Niech S;(ij), S-(i,j) oznaczają odpowiednio elementy tych obrazów 
(piksele) o współrzędnych i,j. Współczynnik korelacji macierzy S7, $2 można zatem zdefinio- 


wać w następujący sposób: 


Me 
Mz 


II 
— 
L. 

II 
= 


(Sili, j) =S) (Sali, J)—(52)) 
= , (3.9) 


A 
Ra 
Ko 


(S,(1,5)-(859)-2, 2 (8,0, J)-48)) 


i=1 j=l 


e! 
Mz 
Me 


II 
= 
p 

Il 
pa 


gdzie symbol <> oznacza uśrednianie według następującego wzoru: 


(S= 2 251). (3.10) 
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Na bazie współczynnika korelacji zbudowany został algorytm wyszukujący lokalne prze- 
mieszczenia. Aby oszacować przemieszczenie wybranego fragmentu z pierwszej fotografii 
wyszukiwany jest najbardziej podobny do niego obraz z drugiej fotografii. W tym celu musi 
on zostać porównany ze zbiorem obrazów odpowiadającym wszystkim potencjalnym prze- 
mieszczeniom. Ponieważ w badanych przepływach możliwe przemieszczenia są ograniczone 
do pewnego zakresu (zależy to od skali fotografii, odstępu czasowego pomiędzy zdjęciami), 
porównaniu podlegają tylko te fragmenty, które odpowiadają przemieszczeniom w danym 


zakresie. 


Należy pamiętać o tym, że obraz cyfrowy nie jest ciągły. Odpowiada on skończonej macie- 
rzy. Można wyodrębnić z niego skończoną liczbę podmacierzy o zadanym rozmiarze. Zatem 
zbiór obrazów odpowiadających wszystkim potencjalnym przemieszczeniom jest ograniczo- 


ny. 


3.2.1 Dodatkowe algorytmy zastosowane w programie PIV-Kor 

Oprócz opisanego powyżej algorytmu wyszukiwania lokalnych przemieszczeń, w progra- 
mie PIV zastosowano szereg pomocniczych algorytmów, które usprawniają jego pracę 
i naprawiają błędnie wyznaczone przemieszczenia. Oryginalnym wkładem autora w metody- 
kę analizy PIV jest algorytm redukujący główny ruch. Dzięki temu program PIV-Kor na wstę- 
pie potrafi zmierzyć główne przemieszczenie (w przypadku eksperymentu w komorze jest to 
ruch opadający struktury chmurowej) i wyeliminować je poprzez odpowiednią modyfikację 
analizowanych obrazów przepływu. Tak spreparowane obrazy są dalej przekazane algorytmo- 
wi szukającemu lokalnych przemieszczeń, ograniczając się tylko na znajdowaniu fluktuacji 


przemieszczeń. 


Elementem poprawiającym dokładność analizy pola prędkości jest również procedura 
iteracyjnej deformacji obrazów. Idea tego postępowania, znana w literaturze analizy obrazów 


[31], została zaimplementowana przez autora w programie PIV-Kor. 


Redukcja głównego ruchu 

W opisywanych w niniejszej pracy eksperymentach ruch główny spowodowany opada- 
niem struktur chmurowych był dużo większy od fluktuacji prędkości. Skutkuje to bezpośred- 
nio tym, że para fotografii wydaje się zawierać ten sam obraz tylko przesunięty. Bezpośrednie 
zastosowanie algorytmu do znajdowania lokalnych przemieszczeń (rozdz. 3.2) wymagałoby 


użycia dużego rozmiaru obszaru przeszukiwania. To skutkowałoby wydłużeniem czasu obli- 
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czeń i zmniejszało dokładność analizy. Tych problemów udało się uniknąć przez specjalną 
konstrukcję programu, który przed przystąpieniem do znajdowania lokalnych przemieszczeń 
eliminuje ruch średni z pary fotografii. Przesunięcie średnie pomiędzy dwiema fotografiami 
jest oszacowane przez zastosowanie podobnego algorytmu jak przedstawiony w par. 9.2. Po 
znalezieniu przemieszczenia obie fotografie są na siebie nasuwane tak, aby widoczne na nich 
struktury jak najlepiej się pokrywały, a „wystające” obszary na brzegach fotografii są odrzu- 
cane. Po tej operacji otrzymuje się fotografie mniejszych rozmiarów niż oryginalne. Odcięcie 
niepokrywających się fragmentów nie pociąga za sobą utraty cennych informacji. Odcinane 
są obszary, które zawierają kropelki występujące tylko na jednej z dwóch fotografii, więc 
zupełnie nieprzydatne do poszukiwania przemieszczeń, a nawet szkodliwe, gdyż mogą być 


przyczyną błędów przy zastosowaniu algorytmu szukającego lokalnych przemieszczeń. 


Wyznaczona przez algorytm wartość przesunięcia jest zapamiętywana i dodawana pod 
koniec obliczeń do uzyskanego pola prędkości. W ten sposób uzyskiwane są całkowite pola 


prędkości w analizowanym przekroju przepływu. 


lteracyjna metoda deformacyjna 


Kolejnym elementem usprawniającym działanie zbudowanego programu jest zniekształca- 
nie obrazu zgodnie z polem przemieszczeń. Pozwala to na znaczne zwiększenie dokładności 


analizy PIV. 


W przypadku przepływów o dużych gradientach prędkości założenie o tym że próbka foto- 
grafii nie zmienia się może być niedostatecznie spełnione. Nie dotyczy to najczęściej całej 
powierzchni pola prędkości, ale pewnych jego części takich jak obszary Ścinania, gdzie 
występują duże gradienty prędkości. Dobre efekty w takich przypadkach może dać dodatko- 
we zniekształcanie obrazów dopasowujące je tak aby zminimalizować występujące różnice 
położeń cząstek i uwzględnienie transformacji tego zniekształcenia przy obliczaniu ostatecz- 


nego wektora przemieszczeń. 


Niech będzie dany wektor przemieszczenia d zaczepiony w punkcie (x,y). Wyznacza 
on przemieszczenie otoczenia tego punktu. Z drugiej strony wektor —d wyznacza przemiesz- 
czenie punktu d+(x,y) z drugiej fotografii po odwróceniu czasu (patrz rys. 3.3). Po znie- 
kształceniu całej drugiej fotografii zgodnie z odwróconym polem przemieszczeń otrzymaliby- 


śmy pierwszą fotografię. 
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d+(x,y) 


(x,y) 


Rys. 3.3: Ilustracja idei metody deformacyjnej. Wektor d jest wektorem przemieszczenia cząsteczki, 
która w chwili początkowej znajdowała się w punkcie (x,y). W chwili wykonania drugiej fotografii 
cząsteczka ta znajdowała się w punkcie d+(x,y). Jeśli przesunie się cząsteczkę z drugiej fotografii o 
wektor -d to znajdzie się ona w swoim początkowym położeniu, czyli tam gdzie znajdowała się na 
pierwszej fotografii. Jeżeli posiadamy kompletne pole przemieszczenia to stosując tę operację do 
wszystkich punktów drugiej fotografii otrzymamy pierwszą fotografię. 


W wyniku zastosowania algorytmu wyszukiwania lokalnych przemieszczeń otrzymujemy 
pole, które z pewną dokładnością odtwarza rzeczywiste pole przemieszczeń. W szczególności 
duże gradienty prędkości zostają wygładzone. Zniekształcenie drugiej fotografii zgodnie 
z tym polem spowoduje, że nie dostaniemy w rezultacie pierwszej fotografię, ale fotografię, 
która będzie do niej bardziej podobna. Operacja taka powoduje, że obraz tej samej cząsteczki 
z drugiej fotografii zbliża się do obrazu tej samej cząsteczki z pierwszej fotografii. Korzyść 
ztego jest taka, że para składająca się z oryginalnej pierwszej fotografii i zniekształconej 
drugiej służy do oszacowania poprawki do pola przemieszczeń. To pozwala na zwiększenie 
dokładności pomiaru. Poprawione pole prędkości służy ponownie do deformacji drugiej foto- 
grafii i pozwala na znalezienie kolejnej poprawki, która coraz lepiej przybliża pole przemiesz- 
czeń. Ponieważ w kolejnych iteracjach tego procesu obrazy tych samych cząsteczek z pierw- 
szej 1 drugiej fotografii zbliżają się do siebie, do obliczeń używa się coraz mniejszych obsza- 
rów przeszukiwania i mniejszych rozmiarów próbki co jest dodatkowym czynnikiem popra- 


wiającym dokładność pomiarów znacznie poniżej rozdzielczości obrazu wynoszącej | piksel. 


W praktyce jest to realizowane w ten sposób, że pikselowi ze zdeformowanej drugiej foto- 


grafii przypisany jest następujący wektor: 
P;,=(24,55,B”), (3.11) 


Dla wygody można się posługiwać tymi wielkościami niezależnie jako funkcjami wektora 


położenia (x, y) na płaszczyźnie obrazu: 
EJ b B”(x,y). (3.12) 


n jest indeksem numerującym iteracje, x,y to współrzędne danego piksela, xo,vo to współ- 


rzedne danego piksela przed deformacją, B” to jasność danego piksela (opisuje zniekształconą 


44 3 Anemometria obrazowa - PIV 


fotografię). Deformacja stosowana w algorytmie zniekształcającym jest opisana następujący- 


mi przekształceniami: 


B'* (x,y)=B"(x+d;,y+d)), (3.13) 
x (x,y)=xg(x+di,y+d)), (3.14) 
yo (x,y)=yo(x+dz,y+dy). (3.15) 


Przy czym warunki początkowe dla n=0 odpowiadają sytuacji przed deformacją, czyli B’ to 
po prostu oryginalna druga fotografia, a X)=X,yy=Y . Pola skalarne xo, yo służą zatem do tego 
aby każdy przemieszczany piksel obrazu „pamiętał” swoje oryginalne położenie na fotografii 
przed zastosowaniem kolejnych iteracji zniekształcania. Pole wektorowe d”(x ,y)=(d;,d/) 
opisuje pole przemieszczeń oszacowane dla oryginalnej pierwszej fotografii i drugiej fotogra- 
fii zniekształconej w n-tej iteracji — B”. 

Równania 3.13-3.15 można interpretować jako adwekcję trzech pól skalarnych spowodo- 


wana przepływem, którego chwilowe pole prędkości jest opisywane przez d"(x,y). 


x,(x,y),yy(x,y) są współrzędnymi początkowymi elementu obrazu, który w n-tym 
kroku wskutek adwekcji spowodowanej deformacjami znalazł się w punkcie (x,y). Służą 
one do tego, aby każdy piksel z deformowanej drugiej fotografii pamiętał o swoim położeniu 
w kroku n=0, czyli całkowite przemieszczenie tego elementu, oszacowane w n-tym kroku, 
wynosi: 


D” (x, y)=(x07x, yo 7y). (3.16) 


Równania 3.13-3.15 można więc zapisać w postaci: 


B'* (x, y)=B°(x+D}, y+ D'), (3.17) 
Xo (x,y)=x(x+D;,y+D)), (3.18) 
yo (x, y)=yo(x+ Di, y+ D}). (3.19) 


Deformowanie oryginalnych pól jest bardziej korzystne niż deformowanie cząstkowe w kolej- 
nych krokach, ponieważ to drugie niesie za sobą niebezpieczeństwo pojawienia się błędów 


i rozmywania obrazu na skutek interpolacji stosowanej przy zniekształceniach. 


W przygotowanym oprogramowaniu zniekształcanie obrazów jest stosowane pomiędzy 
kolejnymi zmniejszeniami analizowanego fragmentu obrazu. Działający algorytm wraz 
z zaimplementowaniem zniekształcania wygląda następująco: Na początku wstępnie zostaje 


oszacowane pole przemieszczeń, które jest następnie traktowane jako pole deformacji. Zgod- 
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nie z nim deformowana jest druga fotografia, która po tej operacji staje się bardziej podobna 
do pierwszej fotografii. W kolejnym kroku wielkość próbki i obszaru przeszukiwań jest 
zmniejszana i znajdowane są poprawki do pola przemieszczeń. Po wielu takich krokach 
deformowana pierwsza fotografia staje się niemal identyczna z drugą, a pole przemieszczeń 


oszacowane jest z dużą dokładnością i wiernością szczegółów. 


Algorytm ten przetestowano na parach obrazów, które składały się z oryginału i jego 
zdeformowanej kopii otrzymanej przez deformację oryginału zgodnie ze sztucznie wygenero- 
wanym polem przemieszczeń (patrz rozdz. 9.3). Zadaniem programu było znalezienie pola 
deformacji, które posłużyło do wygenerowania takiej pary obrazów. Testy te wykazały, że 
średni moduł różnicy pomiędzy oryginalnym polem przemieszczeń, a polem oszacowanym 
przez program wynosił 0,3 piksela a maksymalnie około 1 piksela, gdzie średnie przemiesz- 


czenie wynosiło 7,6 piksela. 


Filtr medianowy 


W praktyce w pomiarach PIV typowa jest sytuacja przedstawiona schematycznie na 
rys. 3.4, gdzie w uzyskanym polu prędkości pojawiają się błędne wektory powodujące 
powstanie nieciągłości pola. Wspomniana wcześniej deformacja obrazu zastosowana zgodnie 
z takim polem prowadziłaby do nieciągłych deformacji i propagacji błędu w kolejnych itera- 
cjach. Problem ten jest znany od początku techniki PIV i aby go uniknąć stosowane są filtry 
medianowe ([30], [32]), które pozwalają na eliminację błędnych wektorów. 


Filtr sprawdza każdy wektor otrzymanego pola prędkości pod kątem jego koherencji 
z wektorami z otoczenia. Załóżmy, że do sprawdzenia koherencji wektora v. branych jest 
p wektorów z jego najbliższego otoczenia. Łącznie z wektorem v. otrzymujemy zbiór o licz- 
bie elementów p+1. Niech elementy tego zbioru będą numerowane indeksem i = 1, 2, ..., p+1 
Zdefiniujmy średnią różnicę pomiędzy wektorem vj; i resztą wektorów ze zdefiniowanego 
zbioru: 

DA | 
A= ; (3.20) 
p 

Dalej możemy znaleźć wektor medianowy, którego indeks powinien spełniać następującą 

zależność: 


med =arg min 4, (3.21) 


l 
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Rys. 3.4: Fragment zmierzonego pola prędkości z jednym błędnym wektorem. 

Jeśli wszystkie wektory z rozważanego zbioru zaczepić w jednym punkcie, to końce 
wektorów utworzą chmurę punktów, której rozciągłość jest charakteryzowana przez wartość 
Amea. Wektor Vmed jest wektorem który znajduje się najbliżej środka chmury. Amea jest więc 
średnią odległością od środka chmury. Aby przetestować koherencję danego wektora 
v. sprawdzany jest następujący warunek: 


NV zad SA nea (3.22) 


Spełnienie tego warunku oznacza, że testowany wektor znajduje się w obrębie klastra 


wektorów, jest więc z nimi koherentny. 
W przeciwnym wypadku, gdy spełniony jest dopełniający warunek: 
NWA: (3.23) 


to wektor v. nie jest koherentny z wektorami otoczenia, a więc z dużym prawdopodobień- 
stwem jest to wektor błędny. Aby naprawić ten błąd, w miejsce wektora v. jest wstawiany 


wektor Vmed. 


3.3 Stereo PIV 


Pole prędkości charakteryzują zawsze trzy jego składowe (U, V, W). Głównym ogranicze- 
niem metody 2D-PIV, która została opisana wcześniej jest to ,że tak naprawdę otrzymujemy 
rzut pola prędkości płynu w płaszczyźnie świetlnej na powierzchnię fotografii. Tracona jest 
przez to informacja o składowej pola prędkości prostopadłej do tej powierzchni. Technika 


Stereo-PIV [33] powstała, aby przezwyciężyć to ograniczenie. 
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Dzięki zastosowaniu zestawu dwóch zsynchronizowanych, identycznych kamer skierowa- 


nych na badany obszar pod różnymi kątami, system stereo może „widzieć” przepływ stereo- 


skopowo. 
Poruszającysię _ 7 
obiekt w dwóch 255 Arnie 
chwilach czasu es 4 \ 
\ s 
1 , 
w \ 
Płaszczyzna Płaszczyzna 
fotografii zlewej fotografii z 
kamery prawej kamery 


Rys. 3.5: Ogólna idea pomiaru stereo PIV. Poruszający się obiekt jest rzutowany przez układy 
optyczne kamer na płaszczyzny fotografii. Ponieważ kamery ustawione są pod różnymi kątami, 
rejestrują one różne rzuty wektora przemieszczenia. Gdy znana jest geometria układu pomiarowego, 
możliwe jest odtworzenie rzeczywistego wektora przemieszczenia z odpowiedniego złożenia dwóch 
rzutów tego wektora. Proporcje na rysunku zostały zniekształcone dla lepszego zobrazowania 
geometrii. 


W technice Stereo PIV rejestruje się, podobnie jak w klasycznej technice, pary fotografii, 
przy czym każda z kamer rejestruje oddzielną parę symultanicznie. Para z każdej kamery jest 
przetwarzana i uzyskiwane są dwa dwuwymiarowe pola prędkości, które następnie składane 
są w jedno pole prędkości zawierające wszystkie trzy składowe. Informacje o geometrii 
pomiaru, potrzebne do złożenia dwóch pól, pobierane są w procesie kalibracji układu, 
w którym wykonuje się fotografie specjalnego wzorca (rys 3.6) umieszczanego w płaszczyź- 
nie świetlnej. Analiza fotografii wzorca pozwala na znalezienie przekształceń z układu odnie- 
sienia fotografii do rzeczywistego układu odniesienia badanego przepływu. Ogólny schemat 


działania techniki Stereo-PIV przedstawiony jest na rys. 3.9. 


Używany w omawianych badaniach techniką Stereo-PIV układ laboratoryjny był podobny 
do układu stosowanego w klasycznej metodzie 2D-PIV. Podstawową różnicą jest zastosowa- 
nie dwóch kamer ustawionych pod kątem 90° względem siebie a pod kątem 45° względem 
płaszczyzny świetlnej. Tak duże kąty obserwacji stwarzają jednak problemy z zapewnieniem 
ostrości obrazu w całej płaszczyźnie. Konieczne jest zastosowanie tzw. adaptera Sheimpfluga 


umożliwiającego skręcenie osi optycznej obiektywu względem osi optycznej kamery. Przy 
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ustawieniu odpowiedniego kąta skręcenia możliwe jest otrzymanie fotografii o dobrej ostrości 
w każdym punkcie fotografii. W doświadczeniach używane były adaptery Sheimpfluga firmy 
ILA. 


Do analizowania otrzymanych w czasie eksperymentów fotografii użyto oprogramowania 
Vid-PIV firmy ILA. Przykładowe pole prędkości struktury chmurowej prezentowane jest na 
rys. 3.10. 


Zasadniczym celem wykorzystania techniki Stereo-PIV w obecnych badaniach była wery- 


fikacja izotropowości fluktuacji turbulentnych pola prędkości w płaszczyźnie poziomej. 


Rys. 3.6: Wzorzec używany do kalibracji układu stereo PIV. Widoczne są białe krzyżyki, których 
położenie jest wykrywane przez program kalibrujący. Dzięki temu, że markery są ułożone na 
regularnej siatce o znanych odstępach pomiędzy oczkami siatki, możliwe jest oszacowanie 
przekształcenia z układu odniesienia fotografii do współrzędnych rzeczywistych. 


3.3 Stereo PIV 


SS 
Rys. 3.7: Układ doświadczalny z systemem Stereo-PIV. Na pierwszym planie widoczny laser 
impulsowy. Wiązka laserowa jest kierowana do szczeliny w komorze wykorzystując tzw. ramię 
Świetlne. Płaszczyzna świetlna tworzona przez zastosowanie układu optycznego z soczewką 
cylindryczną oświetla wnętrze komory, w której zachodzi badany proces mieszania. Przed komorą 
widoczne są kamery w ustawieniu stereoskopowym oraz komputer sterujący systemem i archiwizujący 
dane. 


Rys. 3.8: Widok kamer z adapterami Sheimpfluga w ustawieniu stereoskopowym rejestrujących 
struktury tworzone przez opadającą mgłę w komorze chmurowej. 
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Rys. 3.9: Schemat działania układu Stereo-PIV. 


3.3 Stereo PIV 


Rys. 3.10: Widok perspektywiczny trzech składowych wektorów prędkości zmierzonych dla struktury 
chmurowej w komorze. Barwa od czerwonej do niebieskiej odpowiada wartości składowej prędkości 
prostopadłej do płaszczyzny świetlnej. 


http://rcin.org.pl 
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4 Eksperyment 


4.1 Opis układu eksperymentalnego 


Używany w pracy układ eksperymentalny jest 
rozbudowaną wersją układu używanego wcześniej 
do badań laboratoryjnych zawiesiny chmurowej 
([18], [35], [36]). Głównym elementem układu 
laboratoryjnego jest komora o wymiarach 
lm x Im x 1,8m. Przednia ściana komory wyko- 
nana została ze szkła tak, żeby można było 


swobodnie śledzić zachodzące wewnątrz procesy. 


W bocznej ścianie wykonana została szczelina 
Rys. 4.1: Uproszczony schemat układu 


umożliwiająca wprowadzenie do wnętrza płasz- eksperymentalnego: 1 - mata komora 
"= ‘ a z generatorem kropelek; 2 - główna komora; 
czyzny świetlnej (patrz rys. 4.2). Jedna ze ścian 3 - płaszczyzna świetlna; 4 - laser; 5 - powietrze 
pełni również rolę drzwi, dzięki czemu możliwe CANTONE Miej ZZ 7 REM 
jest otwarcie komory. Na górze głównej komory znajduje się mniejsze pudło o wymiarach 
0,3m x 0,6m x 0,5m. Jest ono napełniane kropelkami wody produkowanymi przez ultradźwię- 
kowy generator. Napełnianie trwa 1,5 minuty, tak żeby uzyskać maksymalną ilość kropelek 
wody w danej objętości powietrza. Oszacowana wodność bezwzględna! (patrz paragraf 4.3) 


powietrza chmurowego używanego w eksperymentach wynosi 24 g/m’. 


Pomiędzy małą komorą a główną komorą znajduje się okrągły otwór o średnicy 10 cm. 
Przegroda zamykająca go pozwala na ustalenie początkowych warunków eksperymentu, 
w których obie masy powietrza (powietrze chmurowe na górze i powietrze czyste na dole) są 
odseparowane. Po otwarciu przegrody powietrze chmurowe opada w dół wskutek zmiany 
gęstości spowodowanej przez parowanie kropelek. Dodatkowo przepływ wymuszany był 
przez pracujący ciągle generator kropelek. 


System termopar i czujników wilgotności umieszczony wewnątrz komory umożliwiał 


monitorowanie warunków termodynamicznych w trakcie eksperymentu. 


I  Wodność bezwzględna — masa kropel wody z których składają się chmury w jednostce objętości wyrażana 
zazwyczaj w jednostkach: [g/m*] [34]. 
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Rys. 4.2: Schemat umocowania kamery i układu formującego wiązkę do prowadnicy. Rzuty płaskie 
układu eksperymentalnego: 


I — rzut z kierunku szczeliny przez którą do komory kierowana jest płaszczyzna Świetlna, 

II — rzut z kierunku przedniej szyby, 

III — rzut z góry. 

1 — główna komora, 2- mała komora, 3 — generator kropelek, 4 — prowadnica, 5 — kamera, 6 — układ 
formujacy wiązkę połączony przegubowym systemem optycznym z laserem (ramie Świetlne — por. rys. 
4.3), 7 — szczelina w głównej komorze umożliwiająca wprowadzenie płaszczyzny świetlnej do wnętrza. 


4.2 Pomiar pól prędkości 


Mieszanie turbulentne obu mas powietrza jest wizualizowane przez zastosowanie światła 
laserowego. W eksperymentach użyto tandemu dwóch laserów impulsowych Nd:Yag (długość 


fali 532nm) o energii 36mJ na puls. Przy zastosowaniu układu optycznego w skład którego 
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wchodzi soczewka cylindryczna możliwe jest uformowanie promienia świetlnego w smugę 
świetlną o grubości około lmm. Smuga ta kierowana jest do wnętrza komory. Przy wyłącze- 
niu innych źródeł światła widoczne stają się tylko te kropelki, które znajdują się w smudze 
świetlnej. Tak zwizualizowane obrazy badanego procesu są rejestrowane przez kamerę CCD 
(lub dwie kamery) i archiwizowane przy użyciu komputera. Do pomiarów pól prędkości 
kropelek użyte zostały kamery PCO Sensi-Cam o rozdzielczości 1280x1024. Kamery te 
sprzężone z tandemem laserowym pozwalają na rejestrację par fotografii, w których odstęp 
czasowy pomiędzy poszczególnymi fotografiami w parze może być regulowany. W przypad- 


ku eksperymentu w komorze chmurowej odstęp ten wynosił 5ms. 
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Przetwomik 


an. - cyfr. 

czujników 
ramie prowadnica wilgotności i 
świetlne temperatury 


Komputer PC 


kamera 


układ formujący 
wiązkę 


Synchronizator 


Rys. 4.3: Schemat urządzeń podłączonych do układu eksperymentalnego. 


Aby usprawnić wykonanie eksperymentów w różnej odległości od wlotu do komory układ 
formujący wiązkę światła i kamerę umieszczono na ramieniu przymocowanym do pionowej 
prowadnicy, tak aby znajdowały się na tym samym poziomie (patrz rys. 4.2). Układ formują- 
cy wiązkę połączony był z laserem przy użyciu ramienia świetlnego co pozwalało na swobod- 
ne zmienianie położenia płaszczyzny świetlnej bez konieczności przemieszczania lasera, przy 


czym płaszczyzna świetlna nie zmieniała swych właściwości. Prowadnica wyposażona 
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w silnik krokowy sterowana była przez komputer. Dzięki temu można było przemieszczać 


w pionie układ rejestrujący bez zmiany położenia kamery względem płaszczyzny świetlnej. 


Ponadto komora została wyposażona w zestaw czujników wilgotności i temperatury sprzę- 
żonych z komputerem. Dzięki temu możliwe było monitorowanie i rejestrowanie temperatury 


i wilgotności wewnątrz komory. 


4.3 Pomiar wodności bezwzględnej 


Istotnym parametrem powietrza chmurowego jest jego wodność. Aby określić jaka ilość 
ciekłej wody zawarta jest w powietrzu chmurowym używanym w eksperymentach wykonano 
pomiary wodności bezwzględnej czyli stosunek masy wody kropelkowej do objętości powie- 


trza w której się znajduje. 


Pomiary polegały na zassaniu z komory chmurowej 10 litrów powietrza przez gumowy 
wąż. Na drodze zasysanego powietrza umieszczono filtr z waty na którym zbierały się kropel- 
ki wody. Obliczając przy użyciu wagi analitycznej (dokładność 0,1mg) różnicę masy suchego 
i mokrego filtru można było określić masę wody kropelkowej znajdującej się w 10 litrach 
powietrza. 
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Rys. 4.4: Pomiary zawartości wody dla różnych wartości masy waty m użytej w pomiarach. Wodność 
wyznaczono dla granicy m0. 
Dla tych samych warunków przeprowadzono szereg pomiarów, których wyniki przedsta- 
wione są na rys. 4.4. Masa filtru była wybierana przypadkowo, dzięki temu można było zaob- 


serwować jak wpływa ona na pomiar wartości wodności. 
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Na podstawie wykresu stwierdzono malejącą liniowo zależność masy zasysanej wody do 
masy filtru użytego w pomiarach. Metodą najmniejszych kwadratów wyznaczono prostą, 
która przecina oś Y z wartością 23,6 g/m. Zależność tę zinterpretowano jako wpływ oporu 
waty na wydajność zasysania. Większa masa waty osłabia ciąg powietrza i powoduje, że 
mniej kropelek jest zasysanych. Przecięcie osi Y można interpretować jako przypadek ideal- 
nego, bezmasowego filtru, który nie stawia żadnych oporów a jednocześnie zbiera kropelki. 
Dla tego wyidealizowanego przypadku pomiar powinien być najbardziej zbliżony do rzeczy- 
wistej wartości. Korzystając z tych założeń i ekstrapolując znalezioną zależność liniową 


możemy oszacować dla danego przypadku wartość wodności na 24+2 g/m’. 


4.4 Widmo kropel 


Krople wody powietrza chmurowego charakteryzowane są przez widmo kropel, czyli 
rozkład statystyczny ich średnic lub promieni. Dla określenia tego widma, kropelki zbierano 
na niewielkiej warstwie oleju rozciągniętej na szkiełku laboratoryjnym. Tak uwięzione 
kropelki fotografowane były w powiększeniu przy użyciu mikroskopu Nikon ECLIPSE 50i 
i kamery PCO 1200hs. Zakładano, że krople w oleju nie parują i nie ulegają zniekształceniu. 


Otrzymane obrazy kropel były następnie analizowane przez własny program napisany 
specjalnie na potrzeby tego eksperymentu. Typowe obrazy mikroskopowe kropel przedstawia- 
ją jasne dyski z ciemną otoczką (por rys. 4.5). Stworzony program lokalizował na fotografii 
ciemne krążki, które odpowiadały obrazowi kropelek. Następnie tworzona była funkcja, która 
była zależnością promienia wodzącego (w układzie zaczepionym w centrum kropelki) od 
unormowanej całki z jasności pikseli po pierścieniu odpowiadającym temu promieniowi 
(unormowanej - tzn. wartość całki podzielona była przez powierzchnię pierścienia). Zlokali- 
zowanie minimum tej funkcji pozwalało na oszacowanie promienia ciemnej obwódki obrazu 
kropelki. Program analizował automatycznie cały zbiór fotografii kropelek wykonanych 
w szeregu eksperymentów. Tym sposobem stworzono próbkę pomiarów o liczebności powy- 


żej 3000. 
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Rys. 4.5: Przykład fotografii kropelek z naniesioną skalą. Obraz wykonany techniką światła 
przechodzącego przy użyciu mikroskopu Nikon ECLIPSE 50i z obiektywem 20x. 
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Rys. 4.6: Przykład działania napisanego na potrzeby niniejszej pracy programu „droplet_analizer". 
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Rys. 4.7: Typowe widmo masowego udziału różnej wielkości kropelek chmurowych występujących w 
doświadczeniach w komorze chmurowej 
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Znając masę kropelek wody w objętości powietrza, na podstawie widma kropelek można 
oszacować liczbę kropelek wypełniającą tę objętość. Na podstawie danych eksperymental- 
nych oszacowano, że w 1m? znajduje się ok 5,6*10' kropelek. Rys. 4.8 przedstawia liczbę 


kropelek w 1m? w zależności od promienia kropelek. 
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Rys. 4.8: Rozkład liczby kropelek w Im’ powietrza w funkcji ich wielkości. 
Można zauważyć, że zarówno masowo, jak i liczbowo średni wymiar kropelek wynosi 


około 4-5 um. 


Niezależny pomiar widma kropelek generowanych na wylocie nawilżacza ultradźwięko- 


wego wykonane aparaturą HELOS firmy Sympatec potwierdził zmierzone optycznie charak- 


terystyki. 


5 Wyniki pomiarów charakterystyk turbulentnego pola 
prędkości 


W tej części pracy prezentowane będą wyniki pomiarów techniką PIV. Pojedyncza seria 
eksperymentalna składała się z rejestracji 100 par fotografii z częstotliwością 3Hz. Wynikiem 


było uzyskanie 100 pól prędkości w każdej serii pomiarowej. 
Kolejne składowe pola prędkości będą oznaczane literami: 
u — pozioma składowa wektora prędkości leżąca w płaszczyźnie świetlnej, 
v - pozioma składowa wektora prędkości prostopadła do płaszczyzny świetlnej, 
w - pionowa składowa wektora prędkości. 


W eksperymentach z wykorzystaniem techniki 2D-PIV możliwe było zmierzenie tylko 
dwóch składowych u i w, natomiast wszystkie trzy składowe pola prędkości zmierzone zosta- 


ły przy użyciu techniki Stereo-PIV. 


W dalszej części pracy pojawiają się fluktuacje prędkości, które były liczone ze wzoru: 
u, =u,—(u,), (5.1) 


gdzie operacja uśredniania < > oznacza tutaj uśrednianie po pojedynczym polu prędkości. 


Prezentowane w dalszej części pracy wyniki pomiarów charakterystyk turbulencji są wiel- 
kościami uśrednionymi po pojedynczych seriach pomiarowych składających się ze 100 pól 
prędkości. Założono tutaj, że przepływ jest jednorodny w obszarze pomiarowym (około 8cm 
0,lcm x 6cm) i statystycznie stacjonarny w czasie wykonywania pojedynczej serii pomiaro- 
wej (30s). Może być więc charakteryzowany przez wartości uśrednione po obszarze pomiaru i 
po czasie wykonywania jednej serii. To zgrubne założenie pozwala na bezpośrednie porówny- 
wanie pomiędzy różnymi seriami pomiarowymi. 

Wszystkie wyniki pomiarów — obrazy struktur chmurowych i odpowiadające im pola pręd- 
kości zostały zarchiwizowane na nośnikach DVD i są dostępne w Zakładzie Mechaniki 


i Fizyki Płynów IPPT PAN. 
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W trakcie badań z wykorzystaniem techniki PIV wykonano trzy rodzaje eksperymentów: 


Pomiar trzech składowych pola prędkości metodą Stereo-PIV z wykorzystaniem 


kropelek wody 


Eksperymenty wykorzystujące metodę Stereo-PIV zostały przeprowadzone w celu stwier- 
dzenia, czy przepływ w płaszczyźnie horyzontalnej jest izotropowy. W pracy prezentowane są 
wyniki pomiarów z serii składającej się ze 100 par fotografii (rozdz. 5.1). Zaniechano wyko- 
rzystania tej techniki do innych badań ze względu na skomplikowaną procedurę techniczną 
kalibrowania układu, gdzie małe błędy w procesie kalibracji propagują się w wynikach. Meto- 
da ma niewątpliwe zalety w sytuacji gdy chcemy zwizualizować przestrzenną strukturę prze- 
pływu. Jest ona jednak mniej wiarygodna, do obliczeń takich wielkości jak współczynnik 


dyssypacji energii, czy funkcje struktury. 


Pomiar dwóch składowych pola prędkości metodą 2D-PIV z wykorzystaniem krope- 
lek wody 


Większość prezentowanych w tej pracy pomiarów pól prędkości zostało przeprowadzo- 
nych techniką 2D-PIV z wykorzystaniem oprogramowania PIV-Kor). W ramach tych badań 
wykonano 53 serie pomiarowe, każda po 100 pól prędkości z czego trzy serie miały na celu 
określić charakterystyki przepływu bezpośrednio na wlocie powietrza chmurowego do głów- 
nej komory (rozdz. 5.2), a 50 serii w celu zbadania charakterystyk przepływu wewnątrz 
komory (rozdz. 5.3). W badaniu przepływu wewnątrz komory zmieniano odległość obszaru 
pomiarowego od wlotu do głównej komory w zakresie 30cm — 70 cm jak również wilgotność 
względną powietrza w zakresie 20% - 65% (patrz tabela 4). Wartości parametrów przepływu 


wyznaczonych na podstawie tych eksperymentów zebrano w tabeli 8 (patrz rozdz. 9.4). 


Pomiar dwóch składowych pola metodą 2D-PIV z wykorzystaniem nieparujących 
kropelek oleju DEHS 


W pracy prezentowane są również wyniki 4 serii pomiarowych z udziałem zawiesiny 
nieparujących kropelek oleju DEHS (patrz tabela 5). Dla tego systemu zmierzono prędkość 
zawiesiny na wlocie i zbadano pola prędkości dla trzech odległości od wlotu do głównej 


komory (30cm, 50cm, 70cm). 
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Tabela 4: Pozycja pomiaru i wilgotność dla eksperymentów w komorze chmurowej z udziałem 
kropelek wody. Numer eksperymentu zgodnie z tabelą 8 w rozdz. 9.4. 


Nr Pozycja Wilgotność Nr Boscia Wilgotnosé 
eksperymentu, [cm] WC eksperymentu [cm] względna 
[%] oes ewes SM [BEC SM 

I 70 20 26 30 36 

2 70 42 27 30 23 

3 70 49 28 50 23 

4 70 61 29 70 23 

5 60 25 30 70 23 

6 60 20 31 70 23 

7 60 29 32 70 23 

8 60 40 33 70 23 

9 60 62 34 70 23 

10 50 62 35 70 38 

11 50 25 36 50 46 

12 50 32 37 50 23 

13 50 42 38 50 29 

14 50 50 39 50 24 

15 50 61 40 50 33 

16 40 59 41 50 24 

17 40 20 42 30 21 

18 40 32 43 30 21 

19 40 41 44 30 21 
20 40 52 45 30 21 
21 40 46 46 30 32 
22 40 48 47 30 42 
23 30 56 48 30 46 
24 30 65 49 30 51 
25 30 20 50 30 31 
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Tabela 5: Porównanie właściwości wody i oleju DEHS. 


Woda Olej DEHS 
Formuła chemiczna H20 CH s04 
Gęstość 1000 kg/m? 912 g/m? 
Ciśnienie pary nasyconej w 2,3 kPa ape 
temperaturze pokojowej 
Zmierzona w eksperymentach 
typowa masa kropelek w jednostce 3 3 
objętości (wodność właściwa dla 24 gm 8 g/m 
wody) 
Lepkość 895 10° Pas 230 10° Pas 
Napiecie powierzchniowe (25°C) 7,2 10°N/m 3,2 10-2N/m 


5.1 Pomiar trzech składowych pola prędkości metodą Stereo PIV 


Przeprowadzono serię pomiarów trzech składowych pola prędkości z wykorzystaniem 
oprogramowania Vid-PIV firmy ILA. Na podstawie zmierzonych pól prędkości wyznaczono 
histogramy fluktuacji dla każdej ze składowych. Na rysunku 5.1 przedstawiono porównanie 


histogramów trzech składowych pola fluktuacji prędkości. 


Czestosc wystepowania, jednostki umowne 
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Rys. 5.1: Porównanie histogramów trzech składowych fluktuacji pola prędkości (u, v — składowe 
poziome, w — składowa pionowa). 
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Widoczna jest wyraźna różnica rozkładu pomiędzy składowymi poziomymi (u v') a skła- 


dową pionową (w'). Odpowiednie odchylenia standardowe w tym przypadku wynoszą: 
0,=5'10 m/s dla składowej poziomej u, 
o=4:10*m/s dla składowej poziomej v, 
Ow=11-10“m/s dla składowej pionowej w. 


Widoczna jest wyraźna anizotropia kierunku pionowego w stosunku do składowej pozio- 
mej. Przyjęto, że niewielkie różnice w rozkładach horyzontalnych składowych prędkości są 
niefizyczne i wynikają z niedokładności pomiaru prędkości prostopadłych do płaszczyzny 
świetlnej. Z geometrii układu laboratoryjnego wynika, że w płaszczyźnie poziomej przepływ 
nie powinien wyróżniać żadnego kierunku. Różnice mogą też być spowodowane niedoskona- 


łościami w procesie kalibracji układu. 


W dalszej części pracy opierając się na wynikach pomiarów Stereo-PIV i wniosków płyną- 
cych z geometrii układu, będziemy przyjmować, że w płaszczyźnie horyzontalnej przepływ 
jest izotropowy w sensie statystycznym. Ponieważ metoda 2D PIV, która głównie była stoso- 
wana w eksperymentach umożliwia pomiar tylko tych składowych pola prędkości, które leżą 
w płaszczyźnie fotografii, pola prędkości otrzymane w czasie tych badań nie zawierają żadnej 
informacji o składowej v. Założenie izotropii przepływu w płaszczyźnie horyzontalnej impli- 
kuje równość statystycznych właściwości składowych u i v pola prędkości. W dalszej części 
pracy do obliczenia wielkości które wymagają wszystkich trzech składowych pola prędkości 


będziemy korzystać z tego założenia. 


5.2 Pomiar prędkości na wlocie 


Dla określenia prędkości powietrza chmurowego na wlocie do głównej komory przeprowa- 
dzone zostały pomiary techniką 2D-PIV. Wykonano 300 pomiarów pola prędkości na wlocie 


komory. Wykazały one, że średnia prędkość wlotu powietrza wynosi 17 cm/s. 


W bezpośredniej bliskości wlotu przepływ formowany jest przez tarcie powietrza o regu- 
larny brzeg wlotu, powoduje to, że przepływ jest zdominowany przez struktury koherentne. 
Do analizy tych struktur wykorzystano metodę rozkładu w bazie empirycznych funkcji orto- 
gonalnych (POD — proper orthogonal decomposition) opisaną dokładniej w dodatku (patrz 
rozdz. 9.5). Na rys. 5.2 przedstawiono pola wektorowe pierwszych sześciu modów dekompo- 


zycji metodą empirycznych funkcji ortogonalnych. 
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Rys. 5.2: 6 głównych modów otrzymanych analizą POD pola prędkości na wlocie do głównej komory. 
Energia względna poszczególnych modów (stosunek energii zawartej w pierwszym modzie do sumy 
całkowitej energii przepływu) wynosi: a,=0,28122, a:=0,1183, a;=0,084026, a4=0,080371, 
a;=0,056164, as=0,038636. 

Struktury te zapoczątkowują proces przekazywania energii kinetycznej głównego ruchu do 
mniejszych struktur wirowych. W obserwowanym przepływie tarcie powietrza chmurowego 
o ścianki, a następnie mieszanie się z powietrzem wypełniającym komorę stanowi główne 
źródło energii turbulencji. Struktury te nie są tak wyraźne jak struktury koherentne obserwo- 
wane w zastosowaniach technicznych [37]. Jest to spowodowane zapewne zmiennymi warun- 
kami wewnątrz komory, które wpływają również na przepływ na wlocie (pompowanie powie- 


trza do wnętrza komory powoduje powstanie wymuszonej konwekcji). Można zauważyć, że 
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pierwsze 6 modów jest w przybliżeniu symetryczne względem odbicia wzdłuż osi pionowej.' 
Jest to zgodne z założeniem izotropii przepływu w płaszczyźnie poziomej. Główna część 
energii (28%) jest skupiona w pierwszym modzie, który jest związany z oscylacjami prędko- 
ści powietrza chmurowego na wlocie. Mod 3 odpowiada pojawianiu się w przepływie 
gradientu pionowej składowej wektora prędkości wzdłuż kierunku pionowego. Oba mody 
(113) są w przybliżeniu symetryczne względem odbicia wzdłuż osi poziomej. Pozostałe 
mody (2, 4, 5, 6) prezentowane na rys 5.2 odpowiadają prawdopodobnie tworzeniu się dużych 
struktur wirowych. Mody 2 i 4 mają porównywalne energie względne (0,12 dla modu 2 i 0,08 
dla modu 4) i w przybliżeniu symetryczne odbicie modu 2 względem osi poziomej odpowiada 
modowi 4. Podobna sytuacja zachodzi dla modów 5 i 6 o energiach względnych 0,06 i 0,04 


odpowiednio. 


Podsumowując, analiza POD przepływu bezpośrednio na wlocie chmurowego powietrza 
pokazuje fluktuacje prędkości wtłaczanego powietrza i tworzenie się silnych gradientów pręd- 
kości pionowej. Ponadto tworzą się duże struktury wirowe, które zapoczątkowują turbulentną 
kaskadę energii. Struktury te są symetryczne względem odbicia wzdłuż kierunku poziomego 
i pionowego. Wyniki analizy POD są zgodne z założeniem izotropii przepływu w płaszczyź- 


nie poziomej. 


W miarę oddalania się od wlotu chmurowego powietrza, przepływ staje się bardziej 
chaotyczny, w związku z tym trudniej jest doszukać się w nim powtarzalnych struktur wiro- 
wych. Próby analizy POD pól prędkości dla pomiarów wykonanych w większych odległo- 


ściach od wlotu chmurowego powietrza nie są w tej pracy prezentowane. 


Analiza POD była stosowana również przy pracach pomocniczych związanych z porówny- 


waniem oprogramowania PIV [38]. 


5.3 Pomiar właściwości przepływu wewnątrz komory 


W tej części pracy zaprezentowane i omówione zostaną wyniki eksperymentalne otrzyma- 
ne z pomiarów pól prędkości w komorze chmurowej. Pojedynczy eksperyment PIV składał 
się z rejestracji 100 par fotografii (około 30s czasu rejestracji). W pracy poddano analizie pola 
prędkości otrzymane na podstawie 50 eksperymentów przeprowadzonych dla różnych wilgot- 
ności (od 20% do 65%) i różnych odległości od górnej krawędzi komory (od 30cm do 70cm). 
Każdy eksperyment generuje ok. 100 pól prędkości zawierających ok. 10° wektorów. Zapew- 


I Za symetrię względem odbicia będziemy przyjmować również przypadek, gdy po odbiciu mod jest 
identyczny do oryginalnego modu z dokładnością do znaku (zwrotu wektorów). 
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nia to odpowiedni zbiór danych do analizy statystycznej charakterystyk turbulencji w komo- 


rze chmurowej. Zmierzone eksperymentalnie charakterystyki zebrane są w tabeli 8 rozdz. 9.4. 
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Rys. 5.3: Wartości wilgotności i odległości od górnej krawędzi komory dla poszczególnych 


eksperymentów. 


5.3.1 Średnia prędkość opadania 


Oszacowano średnią prędkość opadania przez uśrednienie pionowej składowej wektora 


prędkości po zespole wszystkich pól prędkości uzyskanych w czasie każdego eksperymentu. 
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Rys. 5.4: Średnia prędkość opadania w funkcji wilgotności dla pozycji 30cm od wlotu— gwiazdki i 


70cm od wlotu— kółka. 


Otrzymane wyniki wskazują na spadek prędkości opadania wraz ze wzrostem wilgotności 


powietrza w komorze. Przy mniejszych wilgotnościach kropelki parują intensywniej, co 
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sprzyja większemu ochłodzeniu powietrza z opadającej strugi, przez co staje się ona cięższa. 


Jest to zgodne z wynikami rozważań nad diagramem mieszania (patrz rys. 2.2). 


Wyniki pomiarów średniej prędkości opadania dla dwóch wysokości są pokazane na 


rys. 5.4. Prędkość zmienia się również wraz z odległością od wlotu przyspieszając w trakcie 


opadania. 


5.3.2 Mikroskala Taylora 


Mikroskala Taylora jest pewną skalą charakteryzującą przepływ turbulentny. Interpretuje 
się ją jako odległość, w której funkcja kwadratowa będąca przybliżeniem funkcji korelacyjnej 
pola prędkości przyjmuje wartość zero [39]. Mikroskala Taylora niesie pewną informację 


o kształcie wirów w mniejszych skalach przepływu. 


Na podstawie danych z eksperymentów zostały wyznaczone mikroskale Taylora dla skła- 


dowej u i w pola prędkości, które wyrażają się odpowiednio wzorami: 


= (u ”) 42 — ćw”) 
E Ta 
Ox Oz 


2 2 w 


5.2 
y (5.2) 


gdzie u', w' są fluktuacjami składowych prędkości (por. równanie 5.1). 


mikroskala Taylora dla 
składowej pionowej 


pozycja [cm] wilgotność względna [%] 


Rys. 5.5: Mikroskala Taylora dla składowej pionowej pola prędkości w zależności od wilgotności 
względnej i pozycji pomiaru. Kółka oznaczają dane doświadczalne, a powierzchnia - dopasowany 
metodą najmniejszych kwadratów dwuwymiarowy wielomian drugiego stopnia. Pionowe kreski 
ilustrują różnicę pomiędzy wynikami pomiarów a punktami na dopasowanej powierzchni. Wielomian 
został wykreślony w celu lepszej wizualizacji wyników, nie ma on interpretacji fizycznej. 
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Zmierzone w przeprowadzonych eksperymentach Średnie wartości mikroskali Taylora 
wynoszą: od 4,3mm do 8,9mm dla składowej poziomej u i od 5,4mm do 12mm dla składowej 
pionowej w, natomiast stosunek A„A, w poszczególnych eksperymentach przybierał wartości 
od 1,2 do 1,5. Oznacza to, że we wszystkich badanych przypadkach mikroskala Taylora dla 
składowej pionowej była większa od mikroskali Taylora dla składowej poziomej (patrz tabela 


8, rozdz. 9.4). 


Na rys. 5.5 przedstawiona jest zależność mikroskali Taylora dla pionowej składowej pręd- 
kości od pozycji pomiaru i wilgotności względnej. Widać wyraźnie, że małe wilgotności 
i bliskość wlotu powietrza chmurowego sprzyjają występowaniu małych wartości mikroskali, 
co jest związane prawdopodobnie z rozciągnięciem dolnej granicy przedziału inercyjnego 


turbulencji do mniejszych skal. 


Obliczono liczbę Reynoldsa dla mikroskali Taylora opierając się na następującej defini- 
cji [22]: 


1 

TKE | 10 [z 
Re. = SE 5.3 
ey 7 I ( ) 

gdzie TKE jest energią kinetyczną turbulencji: 
TKE=>|(u*)+(97)+(w?)] = [2(u7)+(w7)], (5.4) 


v - lepkością kinematyczną powietrza, a © - enstrofią. Enstrofia została zgrubnie oszaco- 
wana na podstawie jednej składowej wektora wirowości: 


BY ERC 
N= Hap =F (wy), (5.5) 


gdzie w, jest składową wektora wirowości prostopadłą do płaszczyzny pomiaru — jedyną 
składową, która może być oszacowana na podstawie dwuwymiarowych przekrojów przez 


pole prędkości. 


Stwierdzono, że typowa wartość turbulentnej liczby Reynoldsa występująca w ekspery- 
mentach wynosi 40. Warto zwrócić uwagę, że liczba Reynoldsa zdefiniowana w oparciu o 


skalę głównego przepływu i jego prędkość średnią jest rzędu 10*. 


5.3.3 Współczynnik dyssypacji lepkiej i skala Kołmogorowa 
Ważna charakterystyka przepływu turbulentnego to współczynnik dyssypacji lepkiej wyra- 


żający się wzorem: 
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e=v(|—-|). (5.6) 


Został on oszacowany przy założeniu izotropii przepływu w płaszczyźnie horyzontalnej 
z następującego wyrażenia: 
du] 


T `, (5.7) 


emwy(2 


Zaobserwowana w eksperymentach zależność współczynnika dyssypacji lepkiej od wilgot- 
ności powietrza przedstawiona jest na rys 5.6. Zmierzone w eksperymentach wartości miesz- 
czą się w przedziale od 5*10* m/s? do 2,3*10° m*/s* Wartości te są typowe dla chmur cumu- 
lus przy średniej konwekcji [24]. Potwierdza to poprawność założenia podobieństwa zjawisk 
w mikroskali obserwowanych w laboratorium do rzeczywistych przepływów atmosferycz- 
nych. 

Pomiary współczynnika dyssypacji lepkiej wykonano w różnych odległościach od wlotu 
w komorze chmurowej. Zmiana odległości od wlotu charakteryzuje różne etapy rozwoju prze- 
pływu w komorze związane z coraz większym wpływem zasysanego powietrza suchego 


(zmiana parametru k na rys 2.1). 
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Rys. 5.6: Wyniki pomiarów współczynnika dyssypacji lepkiej dla różnych pozycji pomiaru i 
wilgotności. 
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Dla wszystkich pozycji wyraźnie widoczna jest monotoniczna zależność współczynnika 
dyssypacji lepkiej od wilgotności powietrza, co potwierdza hipotezę o wpływie parowania na 


charakterystyki ruchu turbulentnego. 


Rysunek 5.7 ilustruje zbiorczą zależność współczynnika dyssypacji lepkiej od pozycji 
pomiaru i wilgotności powietrza suchego. Aby lepiej zobrazować monotoniczność rozważanej 
zależności, do danych doświadczalnych dopasowano dwuwymiarowy wielomian trzeciego 


stopnia. 


Oprócz omawianego wcześniej malejącego zachowania współczynnika dyssypacji dla 
ustalonej pozycji eksperymentu, rzucającą się cechą jest niemonotoniczne zachowanie dla 
ustalonej wilgotności. Wyraźnie widać, że dla wilgotności około 20% współczynnik dyssypa- 


cji lepkiej maleje wraz z opadaniem powietrza chmurowego po czym zaczyna rosnąć. 


Dopasowaną powierzchnię należy traktować z dużą rezerwą ze względu na niepewności 
pomiaru współczynnika dyssypacji 1 wilgotności. Jednak dla małych wilgotności, gdzie 
dokładność pomiaru wilgotności jest duża eksperymenty zostały powtórzone kilkakrotnie (dla 
pozycji 30cm, 50cm i 70cm) dając dobrą powtarzalność wyników. Przedstawiona zależność 
współczynnika dyssypacji od pozycji pomiaru przy ustalonej wilgotności około 20% nie 
może być więc przypadkowa i musi wiązać się z naturą badanych procesów a nie z niedosko- 


nałością eksperymentu. 


Duży współczynnik dyssypacji dla pozycji 30cm od wlotu jest spowodowany dwoma 
mechanizmami. Pierwszym jest mechaniczna produkcja turbulencji wskutek rozpadu struktur 
koherentnych tworzących się w trakcie pompowania powietrza przez okrągły otwór. Struga 
powietrza chmurowego, opadając w dół, rozszerza się. Energia kinetyczna turbulencji dyfun- 
duje więc w płaszczyźnie horyzontalnej. Zachodzący przez cały czas proces dyssypacji 
lepkiej powoduje stopniowe zanikanie turbulencji wraz z oddalaniem się od wlotu. Na ten 
proces nakłada się oddziaływanie fluktuacji gęstości spowodowanych parowaniem wody 
kropelkowej, której rozkład przestrzenny jest niejednorodny i poddawany ciągłej ewolucji 
wskutek oddziaływania z turbulentnym polem przepływu. Może to stanowić dodatkowe 
źródło energii turbulencji. Powietrze chmurowe opadając w dół miesza się z czystym powie- 
trzem, co powoduje zmianę proporcji zmieszania. Rozkład gęstości prawdopodobieństwa dla 
proporcji zmieszania w badanym przez nas fragmencie przepływu nie jest znany. Średnia dla 
tego rozkładu powinna być jednak malejącą funkcją pozycji pomiaru, obie masy powietrza 
mieszają się wraz z oddalaniem się od wlotu. Rosnący współczynnik dyssypacji dla dużych 


odległości od wlotu można więc wiązać z monotonicznością diagramu mieszania dla małych 
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wilgotności. Malejące proporcje zmieszania powietrza chmurowego powodują powstawanie 


większych różnic gęstości i produkcję dodatkowej porcji energii turbulencji. 


współczynnik dyssypacji 
lepkiej [m?/s*] 


0.024 — 


0.022 


0.02 ~ 


0.018 ~ 


pozycja [em] ee 
Rys. 5.7: Zależność współczynnika dyssypacji lepkiej od pozycji pomiaru i wilgotności. Kółka 
oznaczają dane doświadczalne, a powierzchnia dopasowany metodą najmniejszych kwadratów 
dwuwymiarowy wielomian trzeciego stopnia. Pionowe kreski ilustrują różnicę pomiędzy wynikami 
pomiarów a punktami na dopasowanej powierzchni. Widoczna jest niemonotoniczność zmian 
współczynnika dyssypacji lepkiej. Poziome osie na rys 5.5 i 5.7 są odwrócone względem siebie ze 


względu na lepszą wizualizację wyników. Należy wziąć to pod uwagę przy porównywaniu obu 
wykresów. 


Wraz ze wzrostem wilgotności zależność współczynnika dyssypacji lepkiej zmienia swój 
charakter i dla wilgotności około 60% staje się monotonicznie malejąca. Zmienia się także 
dynamika tej zależności i różnice pomiędzy pomiarami na różnych wysokościach nie są już 
tak znaczące jak to miało miejsce dla małych wilgotności. Zachowanie to można wiązać 
z „płaskim” przebiegiem diagramu mieszania dla wilgotności zbliżonych do wartości 60%. 
Ilustruje to rys.5.8 który przedstawia porównanie diagramów mieszania i zależności współ- 
czynnika dyssypacji lepkiej od pozycji dla wilgotności 20% i 60%. Chociaż nie widać bezpo- 
średniego związku pomiędzy tymi zależnościami to jednak wyraźnie widoczne jest, że dla 
RH=20% duża zmienność zależności temperatury gęstościowej od proporcji zmieszania 
odpowiada dużej zmienności współczynnika dyssypacji lepkiej w zależności od pozycji 
pomiaru. Z diagramu mieszania wynika, że dla wilgotności równej 60% gęstość nie wykazuje 
zmienności w szerokim zakresie proporcji zmieszania. Wynika z tego, że gęstość zmieszane- 
go powietrza słabo zależy od historii mieszania. To przekłada się na współczynnik dyssypacji 


lepkiej, który również wykazuje słabą zależność od odległości od wlotu powietrza chmurowe- 


go. 
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Rys. 5.8: a) - diagram mieszania dla wilgotności 20% i 60% (T=25C, wodność LWC=24g/m'); b) — 
zależność współczynnika dyssypacji lepkiej od pozycji pomiaru dla wilgotności 20% i 60% uzyskana 
przez dopasowanie wielomianu do danych doświadczalnych. 


Zmierzone w komorze chmurowej pola prędkości pozwalają zbadać izotropowość turbu- 
lencji w badanym przepływie. Dla izotropowej turbulencji spełnione są warunki opisane 
następującym równaniem [39]: 
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au \, dlaj+i. (5.8) 
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Równanie to można przekształcić do następującej postaci: 
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( y \=2, dla j#i. (5.9) 
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Warunek ten jest jednym z często stosowanych testów na lokalną izotropię turbulencji 


przepływu. W przypadku eksperymentu w komorze chmurowej typowa wartość ilorazu: 


I ) 


( dw) 
Oz 


ow 
Ox 


wynosi 2, natomiast wartość ilorazu: 
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ME 


wynosi 1,2. 


Świadczy to o dużym zróżnicowaniu charakteru przepływu w zależności od kierunku 
nawet w małych skalach. Niewątpliwą rolę odgrywa tutaj wpływ pola grawitacyjnego i zwią- 


zane z tym efekty wypornościowe. 
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Omawiana na wstępie skala Kołmogorowa określa najmniejszy rozmiar wirów uczestni- 
czących w ruchu turbulentnym i związana jest z przedziałem dyssypacji. W wirach charakte- 
ryzujących się tą skalą gradienty prędkości są na tyle duże, że istotną rolę w tych skalach 
przepływu odgrywają siły lepkościowe, które wygaszają oddziaływania inercyjne i pochłania- 
ją energię kinetyczną z przepływu. Na podstawie wyników eksperymentalnych oszacowano 
skalę Kołmogorowa korzystając z definicji wyrażającej się następującym wzorem: 


3 1/4 


14 
E 


j= (5.10) 


Typowa dla większości eksperymentów wartość skali Kołmogorowa otrzymana dla lepko- 
ści powietrza i z pomiarów współczynnika dyssypacji lepkiej wynosi 7,6°10*m. Warto zwró- 
cić uwagę, że zmienność współczynnika dyssypacji lepkiej implikuje zmienność skali Kołmo- 
gorowa w funkcji wilgotności i odległości od wlotu chmurowego powietrza (por. rys. 5.7 


15.9): 


skala długości 
Kołmogorowa 
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Rys. 5.9: Skala długości Kołmogorowa w zależności od wilgotności względnej i pozycji pomiaru. 
Kółka oznaczają dane doświadczalne, a powierzchnia - dopasowany metodą najmniejszych 
kwadratów dwuwymiarowy wielomian drugiego stopnia. Pionowe kreski ilustrują różnicę pomiędzy 
wynikami pomiarów a punktami na dopasowanej powierzchni. Poziome osie na rys 5.9 i 5.5 są 
odwrócone względem siebie ze względu na lepszą wizualizację wyników (inny kierunek pozycji 
pomiaru). Należy wziąć to pod uwagę przy porównywaniu obu wykresów. 


Na tys. 5.9 przedstawiono zależność skali Kołmogorowa od parametrów pomiaru — wilgot- 
ności i pozycji. Zależność ta jest jakościowa zgodna z zależnością mikroskali Taylora od tych 


parametrów (por. rys 5.5). Jest to zapewne spowodowane tym, że w obecności kaskady ener- 
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gii, struktura przepływu w najmniejszych skalach (skala Kołmogorowa) zależy od transferu 
energii z większych skal (np. mikroskala Taylora). Skale te są ze sobą bezpośrednio powiąza- 
ne i mniejsza mikroskala Taylora implikuje przedłużenie kaskady energii do mniejszych skal, 


a więc zmniejszenie skali Kołmogorowa. 


5.3.4 Statystyki fluktuacji pola prędkości. 

Podobnie jak przy analizie zmierzonych fluktuacji trzech składowych pola prędkości 
(rozdz. 5.1), wyznaczono histogramy fluktuacji prędkości dla badań wykonanych metodą 2-D 
PIV, mając do dyspozycji składową poziomą i pionową pola prędkości. 

Rysunek poniżej przedstawia histogramy dla dwóch składowych pola prędkości dla jednej 


z serii pomiarowych. 


—— składowa pozioma 
= = = składowa pionowa 


Czestosc wystepowania — jednostki umowne 


u, w [m/s] 
Rys. 5.10: Porównanie rozkładów częstości występowania dla składowej pionowej i poziomej 
fluktuacji prędkości dla serii nr 1. 
Dla badanego rozkładu fluktuacji prędkości policzono następujące wielkości: 
odchylenie standardowe: 


o,=Viu") , o =Viw”); (5.11) 


skośność: 


u 3 0 dw 2 (5.12) 
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spłaszczenie: 


Ay KE BEŻ 2 (5.13) 


Typowe, uśrednione wartości tych wielkości otrzymane w eksperymentach zamieszczone 
są w tabeli 6, natomiast tabela 8 (rozdz. 9.4) przedstawia wartości otrzymane w poszczegól- 


nych eksperymentach. 


Tabela 6. Zestawienie wielkości dla statystyk prędkości fluktuacyjnych — dane z doświadcze- 
nia w komorze chmurowej. 


Odchylenie Skośność Spłaszczenie 
standardowe - © S K 
Składowa pozioma 4,6:10* m/s -0,03 0,0 
Składowa pionowa 5,810” m/s 0,01 -0,2 


Większa wartość odchylenia standardowego dla składowej pionowej związana jest 
z rozciągnięciem struktur przepływu w tym kierunku, świadczy też o tym że większa część 
energii kinetycznej skupiona jest w pionowych fluktuacjach pola prędkości. 

W granicach dokładności pomiarów nie stwierdzono natomiast zależności anizotropii 
turbulencji od wilgotności i pozycji w komorze. Zmierzono stosunek wariancji dla składowej 
pionowej pola prędkości do składowej poziomej pola prędkości ( <w'>/<u?>). Średnia dla 


wszystkich pomiarów wynosiła 1,74 przy średnim odchyleniu standardowym 0,3. 


Otrzymane wartości dla skośności i spłaszczenia wskazują, że rozkłady prędkości turbu- 


lencyjnych nie odbiegają od rozkładu normalnego, dla którego S=0, K=0. 


5.3.5 Funkcje struktury i korelacje 

W analizie przepływów turbulentnych istotną rolę odgrywają funkcje struktury i funkcje 
korelacyjne zwane korelacjami. Przestrzenne korelacje prędkości i funkcje struktury określają 
statystyczną współzależność składowych fluktuacji prędkości mierzonych w różnych punk- 
tach przepływu. 

Funkcje struktury drugiego stopnia opisują średnie z kwadratu różnicy prędkości w dwóch 


punktach oddalonych o pewien wektor w zależności od długości tego wektora. 
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Dla danej składowej pola prędkości można zdefiniować wzdłużne i poprzeczne funkcje 
struktury. Wzdłużna — gdy wektor separujący jest równoległy do danej składowej pola pręd- 
kości, poprzeczna — gdy jest prostopadły. 


Dla otrzymanych pól prędkości wyznaczone zostały następujące funkcje struktury: 


Podłużne funkcje struktury drugiego stopnia: 


SY (1)=([w(x,z+1)-w(x,z)f). (5.14) 
Poprzeczne funkcje struktury drugiego stopnia: 

Si (D=((u(x,2z+1)—u(x,z)P), 

Si(1)=([w(x+, z)—w(x, z)P). (5.15) 


Gdzie / jest długością wektora separującego, a u i w składowymi pola prędkości. 


Na rys. 5.11 i 5.12 można zauważyć, że funkcje struktury zachowują się jak Ý dla małych / 
w okolicach skali Kołmogorowa. Jest to zachowanie typowe dla skal przedziału dyssypacyj- 
nego [40]. Dla 10°m funkcja struktury zmierza do [°° co jest zarazem typowe dla obszaru 
inercyjnego [22]. Widoczna jest anizotropia objawiająca się różnicami w przebiegu poprzecz- 


nych funkcji struktury (rys. 5.12). 


I [m] 


Rys. 5.11: Logarytmiczne wykresy, podłużnych funkcji struktury z serii nr l. Przerywana linia — 
składowa horyzontalna prędkości, ciągła — składowa pionowa. 
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Rys. 5.12: Logarytmiczne wykresy, poprzecznych funkcji struktury z serii nr.l. Przerywana linia — 
składowa horyzontalna prędkości, ciągła — składowa pionowa. 


Rys. 5.13 przedstawia funkcje struktury otrzymane z pomiarów wykonanych przy małych 
wilgotnościach powietrza. Porównane zostały funkcje struktury na różnych wysokościach. 
(30cm, 50cm i 70cm od wylotu). Widoczne są istotne różnice, które świadczą o tym, ze zmie- 
niają się właściwości przestrzenne przepływu. Największe wartości przyjmuje funkcja struk- 
tury dla pozycji 30cm co świadczy o tym, że w tym obszarze pole prędkości charakteryzuje 
się największą zmiennością przestrzenną i wirami o dużej energii co przekłada się na duże 
wartości współczynnika dyssypacji lepkiej. Porównanie podłużnych funkcji struktury dla 
składowej pionowej wskazuje na znacznie większe wartości dla pozycji 30cm od pozostałych 
pozycji pomiarowych. Jest to związane z występowaniem silnych gradientów pionowych 


wskutek dużych przyspieszeń, które występują w tym rejonie przepływu. 


Poprzeczne funkcje struktury są związane z gradientami pola prędkości, które występują 
na przekroju przez struktury wirowe. Dają one zatem informację o przestrzennych własno- 
ściach wirów. Rys.5.13 c) i d) przedstawiają poprzeczne funkcje struktury odpowiednio dla 
składowej poziomej i pionowej. Na rys. c) funkcja struktury dla składowej poziomej dla 
pozycji 70cm przyjmuje większe wartości niż dla pozycji 50cm. Natomiast na rys d) funkcja 
struktury dla składowej pionowej dla pozycji 50cm przyjmuje większe wartości. Możliwą 
interpretacją tego faktu jest to, że w różnych etapach przepływu struktury wirowe podlegają 
deformacjom, które polegają na rozciąganiu bądź zgniataniu wirów wzdłuż kierunku piono- 
wego. Można więc wnioskować, że struktury wirowe są w pozycji 50cm rozciągane, nato- 


miast w pozycji 70cm raczej zgniatane wzdłuż kierunku pionowego. 
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Rys. 5.13: Porównanie funkcji struktury na różnych wysokościach: a) podłużne funkcje struktury dla 
składowej horyzontalnej prędkości, b) podłużne funkcje struktury dla składowej pionowej prędkości, 
c) poprzeczne funkcje struktury dla składowej horyzontalnej prędkości, d) poprzeczne funkcje 
struktury dla składowej pionowej prędkości. Pomiary uśrednione dla serii wykonanych przy małej 
wilgotności (20%-23%). 


Rys. 5.14 przedstawia porównanie podłużnych funkcji struktury dla składowej horyzontal- 
nej i pionowej. Dla wszystkich pozycji pomiarowych widoczne jest, że w całym zakresie skal 
większe wartości przyjmują funkcje dla składowej pionowej prędkości, co najwyraźniej 
widoczne jest na przykładzie pomiarów dla pozycji 30cm. Świadczy to o występowaniu 
dużych gradientów prędkości pionowej wzdłuż kierunku przepływu średniego. Są one zwią- 
zane głównie z geometrią przepływu — ruchem średnim wzdłuż kierunku grawitacji, fluktu- 


acjami siły wyporu wywołanymi głównie procesem mieszania i parowania. 


Analogicznie rys. 5.15 przedstawia porównanie poprzecznych funkcji struktury dla obu 
składowych. Widoczna jest jeszcze większa anizotropia niż to miało miejsce dla wzdłużnych 
funkcji struktury. Świadczy to o występowaniu w przepływie dużych różnie składowej piono- 
wej prędkości wzdłuż kierunku poziomego. Można to wyjaśnić przez rozciąganie struktur 
wirowych w kierunku pionowym spowodowane przez grawitację i siły wyporu. Interesujące 


jest, że efekt ten występuje w całym zakresie wektora separacji /. 
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Rys. 5.14: Porównanie podłużnych funkcji struktury dla składowej poziomej (linia ciągła) i pionowej 
(linia przerywana) pola prędkości: a) — pozycja 30cm, b) — pozycja 50cm, c) — pozycja 70cm. 
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Rys. 5.15: Porównanie poprzecznych funkcji struktury dla składowej poziomej (linia ciągła) i pionowej 
(linia przerywana) pola prędkości: a) — pozycja 30cm, b) — pozycja 50cm, c) — pozycja 70cm. 
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Zaobserwowano również zmianę funkcji struktury wraz z wilgotnością. We wszystkich 
przypadkach dla większych wilgotności funkcje przyjmują mniejsze wartości (patrz 
rys. 5.16). Świadczy to więc o tym, że rosnąca wilgotność tłumi zmienność przestrzenną pola 
prędkości. 
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Rys. 5.16: Porównanie podłużnych funkcji struktury dla różnych wilgotności (pozycja 30cm). 


Kolejnym narzędziem analizy struktur turbulentnych jest funkcja autokorelacji pola pręd- 
kości. Funkcja ta pozwala zbadać statystyczne współzależności prędkości przepływu w 
różnych punktach pola prędkości. Funkcja autokorelacyjna opisuje średnią z iloczynu prędko- 
ści w dwóch punktach oddalonych o pewien wektor w zależności od długości wektora separu- 
jącego oba punkty (Średnia liczona jest dla wszystkich par punktów oddalonych o dany 
wektor). Analogicznie do funkcji struktury, dla danej składowej pola prędkości można zdefi- 
niować wzdłużne i poprzeczne funkcje korelacji. Wzdłużna — gdy wektor separujący jest 


równoległy do danej składowej pola prędkości, poprzeczna — gdy jest prostopadły. 


Dla długości / wektora separacji korelację wzdłużną można zdefiniować w następujący 


sposób: 
R? (1)=(u'(xtl,z)-u'(x,z)), (5.16) 
Ri (1)=(w' (x, z+z)-w'(x,z)). ` 
Analogicznie definiujemy korelacje poprzeczne. 
R;()=lu'(x,z+l)su'(x,z)), (5.17) 


RĘ(D=(w'(x+1,z):w'(x,z)). 
Po unormowaniu przez odpowiednie odchylenia standardowe fluktuacji prędkości otrzy- 


mane funkcje autokorelacji przyjmują wartości jedynie z przedziału [-1:1]. 


Na wykresach korelacji (rys. 5.17 i rys. 5.18) również widoczne są różnice między krzy- 


wymi obliczonymi dla różnych składowych. Szczególnie jest to widoczne dla korelacji 
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podłużnych. Duże wartości korelacji podłużnych dla składowej pionowej mogą być związane 


między innymi z występowaniem gradientu prędkości pionowej wzdłuż tego kierunku. 


Na wykresie korelacji poprzecznych krzywa dla składowej pionowej przyjmuje dla 
1>0,025m wartości ujemne. Jest to związane z występowaniem wyraźnych struktur wirowych 
w przedziale inercyjnym propagujących się w kierunku pionowym — korelacje poprzeczne 
przyjmują wartości minimalne dla długości odpowiadającej rozmiarom dominujących wirów 


w przepływie. 
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Rys. 5.17:Wykresy uśrednionych podłużnych funkcji autokorelacji składowych pola prędkości dla serii 
pomiarowej nr. 1. Przerywana linia — składowa horyzontalna prędkości, ciągła — składowa pionowa. 
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Rys. 5.18: Wykresy uśrednionych poprzecznych funkcji autokorelacji składowych pola prędkości dla 
serii pomiarowej nr I. Przerywana linia — składowa horyzontalna prędkości, ciągła — składowa 
pionowa. 
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5.4 Analiza pola prędkości dla nieparujących kropelek 


Dla zbadania hipotezy o wpływie parowania kropel na struktury turbulentne wykonane 
zostały również eksperymenty, w których wykorzystano nieparujące kropelki syntetycznego 
oleju DEHS (Di-Ethyl-Hexyl-Sebacat, C6Hs004). Gęstość tej substancji wynosi 912 kg/m’. 
Ciśnienie pary nasyconej w warunkach, w których wykonywane były eksperymenty jest 


mniejsze od 1Pa, praktycznie więc jest to substancja nieparująca. 
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Rys. 5.19: Widmo promieni kropelek oleju syntetycznego DEHS używanych w doświadczeniach 
w komorze chmurowej. 


Eksperymenty z olejem przeprowadzono w identyczny sposób jak eksperymenty z wodą, 
ztym, że do produkcji zawiesiny wykorzystywano rozpylacz. Użycie ultradźwiękowego 
nawilżacza, który był używany do eksperymentów z wodą okazało się niemożliwe. Widmo 
wielkości kropelek oleju produkowanych w eksperymentach z olejem DEHS przedstawione 
jest na rys 5.19. Jest ono zbliżone do widma kropel chmurowych choć wyraźnie większy jest 


udział kropel o wymiarach z zakresu poniżej 5 um. 


Mniejsza gęstość użytego oleju, nieco inne rozmiary kropelek, a co wydaje się być czynni- 
kiem najważniejszym — brak parowania i związanego z nim negatywnego wyporu powodo- 
wał, że zawiesina DEHS opadała z o wiele mniejszą prędkością niż jej wodny odpowiednik. 
Dla wygenerowania bardziej intensywnego przepływu w kierunku pionowym podjęto próby 
schłodzenia zawiesiny w mniejszej komorze, jeszcze przed jej wpuszczeniem do głównej 
komory eksperymentalnej. Wskutek mieszania z powietrzem efekt ten w niewielkim stopniu 


zmieniał szybkość opadania obłoku. 


Tabela 7 przedstawia podstawowe charakterystyki fluktuacji pola prędkości zmierzonych 
dla tej zawiesiny. Wyniki te wskazują na to, że przepływ bardzo szybko wyhamowuje wraz 


z oddalaniem się od wlotu. Jest to oczywiste gdyż mieszające się powietrze nie potęguje 
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procesu parowania a przeciwnie, podnosi temperaturę obłoku i zmniejsza koncentrację krope- 
lek co prowadzi do zmniejszenia efektów wypornościowych. Współczynnik dyssypacji 
również monotoniczne zanika z odległością od przegrody łączącej obie komory. To w istotny 
sposób różni się od przypadku z parowaniem kropelek wody. W przypadku przepływu 
z zawiesiną wodną, lokalne ochłodzenie niejednorodnie zmieszanego powietrza powodowało 
powstanie różnic sił wyporu, które z kolei powodowały powstanie w komorze ciągu pionowe- 
go i przyspieszanie strugi opadającego powietrza. Oprócz tego małe skale fluktuacji wyporu 
oddziaływały bezpośrednio z fluktuacjami pola prędkości, co wpływało na energię turbulencji 
i w konsekwencji na współczynnik dyssypacji lepkiej. W przypadku kropelek oleju syntetycz- 
nego współczynnik dyssypacji zanika monotonicznie wraz z zanikaniem prędkości średniej 


przepływu, nie ma źródła podtrzymującego generację turbulencji. 


Tabela 7: Wyniki pomiarów dla eksperymentów z DEHS 


Odległość od Prędkość Współczynnik Stosunek 
wlotu opadania dyssypacji wariancji 
[cm] [m/s] lepkiej [m’/s*] | <w'>/<u> 
30 0,026 0,00020 1,15 
50 0,018 0,00006 1,18 
70 0,009 0,00001 1,17 


<w">/<u">, tabela 7). 


Wyniki pomiarów wskazują również na to, 


że przepływ nie jest izotropowy (patrz 


x10 
5 T T T T T 
ZE S; [m‘/s*] 
45 Il 2/224] | 
S, [m/s] 
4p mi 
3.5 4 
z 3 1 
2 
“E 2.5 4 
© ol ] 
1.5F 4 
1k 4 
0.5 
0 1 i L 1 1 ri 
0 5 10 15 20 25 30 35 40 


| [mm] 


Rys. 5.20: Badanie obłoku oleju DEHS dla pozycji pomiaru 50cm. Porównanie podłużnych funkcji 
struktury dla składowej horyzontalnej (ciągła linia) i pionowej pola prędkości (przerywana linia). 
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Rys. 5.21: Badanie obłoku oleju DEHS dla pozycji pomiaru 50cm. Porównanie poprzecznych funkcji 
struktury dla składowej horyzontalnej (ciągła linia) i pionowej pola prędkości (przerywana linia). 


Porównanie funkcji struktury dla badań zawiesiny DEHS (rys.5.20 i 5.21) wskazuje na to, 
że podobnie jak to miało miejsce dla eksperymentów z wodą (patrz rys. 5.14 5.15), przepływ 


turbulentny z kropelkami oleju charakteryzuje się anizotropią w szerokim przedziale skal. 


Pomimo istotnych różnic między eksperymentami z kropelkami wody i oleju DEHS, doko- 
nano jakościowego porównania wyników z obu realizacji eksperymentu. W obu przypadkach 
widoczna jest anizotropia turbulencji, która wynika z wyróżnionego kierunku przepływu Sred- 
niego i grawitacji. Opadając w komorze, przepływ z wodą jest przyspieszany (patrz rys. 5.4), 
co jest spowodowane parowaniem i zmianą gęstości powietrza, podczas gdy dla przepływu z 
nieparującymi kropelkami oleju obserwujemy gwałtowne wyhamowywanie przepływu (patrz 


tab. 7). 


Współczynnik dyssypacji dla nieparującego oleju (tabela 7) maleje wraz z odległością od 
wylotu (patrz rys. 5.4) co jakościowo odpowiada sytuacji dla eksperymentu z zawiesiną 
wodną przy dużej wilgotności. Można to tłumaczyć tym, że w obu przypadkach zaniedbywal- 


ne są efekty związane z przemianami fazowymi i zawiesina odgrywa rolę pasywną. 


Dla kropelek wody przy małych wilgotnościach współczynnik dyssypacji jest niemonoto- 
niczną funkcją odległości od wlotu co odróżnia go od przepływu z olejem DEHS. To jako- 
Sciowe zróżnicowanie obu przepływów należałoby wiązać z wpływem parowania na strukturę 
przepływu. 

Podsumowując, stwierdzono wpływ parowania kropel na struktury turbulentne przepływu 
w komorze chmurowej objawiające się głównie zmianami współczynnika dyssypacji lepkiej 


izmianami funkcji struktury. Zmierzona w warunkach eksperymentalnych anizotropia 
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w dużej części wynika ze struktury samego przepływu z preferowanym pionowym kierun- 


kiem głównym. 


Niewątpliwie w rzeczywistych chmurach ruch w kierunku pionowym jest generowany siła- 
mi wypornościowymi, a więc można przypuszczać, że przy braku innego wymuszenia jest to 


źródłem silnej anizotropii fluktuacji pól prędkości. 


5.5 Dyskusja powtarzalności eksperymentów i dokładności 
pomiarów 

Przy analizie pomiarów eksperymentalnych narzędziami statystycznymi ważna jest 
kwestia powtarzalności doświadczeń. Można sobie wyobrazić chaotyczny przepływ, który 
w powtarzanym wielokrotnie eksperymencie zmieniałby swój charakter mimo zachowania 
identycznych warunków początkowych i brzegowych. Istnieją więc dwie przyczyny dla 
których tak samo przeprowadzone eksperymenty dawałyby znacząco inne wyniki. Pierwszą 
jest sama natura badanego przez nas zjawiska. Istotnie, gdyby porównywać chwilowe pola 
prędkości w eksperymencie z przepływem chaotycznym, otrzymanie podobnych wyników 
jest bardzo mało prawdopodobne, w tym przypadku interesujące są jednak wartości opisujące 
„makroskopowe” własności przepływu uzyskane przez analizę pewnej próbki pomiarów. 
Pojedynczy eksperyment polega w takim przypadku na zebraniu odpowiedniej próbki i okre- 
śleniu jej właściwości takich jak np. uśredniony współczynnik dyssypacji lepkiej, charaktery- 
styczne skale itp. Problem powtarzalności eksperymentu ogranicza się w takim przypadku do 
powtarzalności własności próbek otrzymanych w wyniku realizacji tego samego eksperymen- 


tu. 


Drugą przyczyną, która wpływa na powtarzalność eksperymentów jest ustalenie warunków 
brzegowych i początkowych. W prezentowanych badaniach istotnym problemem było 
dokładne ustalenie wilgotności. Wzrost wilgotności w komorze był uzyskiwany przez wpusz- 
czenie do komory porcji powietrza chmurowego, która odparowując zwiększała koncentrację 
pary wodnej. Metoda ta obarczona jest niestety wieloma niedoskonałościami. Wilgotność nie 
ustalała się w ten sposób na stałym poziomie, ale fluktuowała z powodu dynamiki procesu 
i nieszczelności komory. Czujniki wilgotności umieszczone były z dala od osi komory, 
w której był obserwowany przepływ, tak aby nie powodować jego zakłóceń. W obecności 
nierównomiernego rozkładu wilgotności powoduje to, że mierzona wartość nie musi odpo- 


wiadać dokładnie tej, która jest założona. Przyjęto zatem, że błąd pomiaru wilgotności wynosi 


5.5 Dyskusja powtarzalności eksperymentów i dokładności pomiarów 89 


5%. Jednak, jeśli pomiar był wykonywany w wilgotności równej wilgotności otoczenia wyżej 
wymienione ograniczenia nie występowały. Po wywietrzeniu komory, jej wnętrze napełniało 
się identycznym powietrzem jak to znajdujące się w otoczeniu. W badanej skali gradienty 
wilgotności powietrza znajdującego się w laboratorium są do pominięcia. Zatem po wywie- 
trzeniu komory przestrzenny rozkład wilgotności staje się jednorodny i pomiar nie zależy od 


miejsca. 


Aby zbadać powtarzalność pomiarów wykonano kilka eksperymentów dla wilgotności 
równej wilgotności otoczenia, to znaczy przed każdym pomiarem komora i laboratorium były 
dokładnie wietrzone. Pomiary wykonano dla dwóch różnych położeń: 30cm i 70cm od wlotu 
powietrza chmurowego. Dla pozycji 30cm wykonano 4 eksperymenty przy wilgotności 
względnej 21%, dla pozycji 70cm — 6 pomiarów przy wilgotności względnej 23%. Z otrzyma- 
nych pól prędkości zmierzono współczynnik dyssypacji lepkiej. Otrzymano średnio s=0,0225 
m*/s* dla odległości od wlotu równej 30cm i e=0,0145 m”/s* dla odległości od wlotu równej 


70cm. Przy czym odchylenia standardowe wynoszą odpowiednio: 0,0006 m/s? i 0,0010 m/s*. 


Małe odchylenia standardowe w porównaniu do zmierzonych wartości współczynnika 
dyssypacji lepkiej świadczą o dobrej powtarzalności pomiarów przy ustalonych warunkach. 
Większe odchylenie standardowe dla pozycji 70cm jest zapewne związane z intermittencją 
spowodowaną nieregularnym ruchem strugi powietrza chmurowego względem osi komory. 


Amplituda tego odchylenia zwiększa się wraz z odległością od wlotu. 


6 Badania numeryczne 


6.1 Oddziaływania pojedynczych kropel ze strukturą wirową 


Dla przeanalizowania procesu parowania kropli na otaczający ją przepływ przeprowadzo- 
no badania numeryczne oddziaływań pojedynczych kropel wody z zadanym polem prędkości. 
Model, który opisywałby oddziaływania kropelek z przepływem i wpływ parowania kropelek 


musi uwzględniać między innymi: 
— powietrze jako mieszaninę suchego powietrza i pary wodnej, 
- wymianę energii wewnętrznej, pędu i pary wodnej pomiędzy kropelką a otoczeniem, 
— parowanie kropelek jako źródło pary wodnej, 


— dyfuzję pary wodnej i lokalne zmiany właściwości powietrza spowodowane zmianą 


koncentracji pary wodnej, 


— parowanie kropelki jako lokalna zmiana energii wewnętrznej (parująca woda pobiera 


z otoczenia energię przemiany fazowej), 

- lokalne zmiany pola temperatury spowodowane parowaniem, adwekcją, przewodnic- 
twem cieplnym, 

- lokalne zmiany gęstości powietrza spowodowane zmianą koncentracji pary wodnej 


i zmianami temperatury. 


Dla uproszczenia, oddziaływanie kropelek na przepływ ograniczymy do oddziaływania punk- 
tów materialnych. Zaniedbujemy w tym przypadku zmiany przepływu związane z geometrią 


i z ruchem cieczy w kropelce. 


Model uwzględnia więc trzy składniki: suche powietrze, para wodna i ciekła woda. 
Potrzebne jest więc rozpisanie równań dla tych składników i określenie ich właściwości mate- 
riałowych. Dla ułatwienia, wielkości dla poszczególnych składników będziemy oznaczać 


następującymi indeksami: 
a — powietrze suche, 


v— para wodna, 
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d — kropelka (ciekła woda), 
f-— otaczający gaz czyli mieszanina suchego powietrza i pary wodnej. 
W obliczeniach dynamiki użyty problem został podzielony na dwie części: 
- przepływ powietrza, 
- trajektorie kropelek. 
Skład powietrza jest charakteryzowany przez ułamki masowe: dla pary wodnej Y, i dla 
czystego powietrza Y,, spełniające następujący warunek: 
TRL, (6.1) 
Niech p; oznacza gęstość mieszaniny powietrza suchego i pary wodnej, wtedy gęstość skład- 
nika j oblicza się z ułamku masowego w następujący sposób: 
P;=Y pep. (6.2) 
Przepływ powietrza liczony był metodą objetości skończonych na podstawie układu 


równań, w którego skład wchodziło prawo zachowania pędu, masy i energii: 


Równanie zachowania pędu 


Zlo) +V) lo v)=-V p+V-T+pyg+F, (6.3) 


Gdzie T' jest tensorem naprężeń lepkich, F siłą zewnętrzną z jaką kropelki oddziaływają 
na przepływ. 
Prawo zachowania masy 


Op 


3, V= Sa * (6.4) 


J 


Sm jest źródłem masy — opisuje produkcję pary wodnej w procesie parowania kropelki. 


Równanie energii 


4, CE) +V (Vlo E+ p))=V|EVT-X hyd ATV) +S, (6.5) 


Gdzie k jest współczynnikiem przewodności cieplnej, J; dyfuzyjnym strumieniem skład- 
nika j: 
J;=-p,DNY,, (6.6) 


gdzie D; jest współczynnikiem dyfuzji składnika j a Y; ułamkiem masowym. 


I T,=Ż u(zl ou +] LS ou.) „gdzie u jest współczynnikiem lepkości dynamicznej. 


Ox, Ox; 573 ij Ox, 
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Pierwsze trzy człony odpowiadają za zmianę energii w wyniku przewodnictwa cieplnego, 
dyfuzji molekularnej i dyssypacji. Ostatni człon to zmiana energii spowodowana występowa- 


niem źródła ciepła. 


Energia całkowita wyraża się wzorem: 


2 


E=h Pp , 6.7 


Całkowita entalpia to średnia ważona entalpii wszystkich składników płynu: 
>, (6.8) 
gdzie h; jest entalpia, a Y; ułamkiem masowym j-tego składnika. 
Entalpia poszczególnych składników wyraża się przez ciepło właściwe przy stałym ciśnie- 
niu: 


b] 6:90. (6.9) 


Równanie transportu pary wodnej 
Ułamek masowy pary wodnej zmienia się lokalnie wskutek konwekcji, dyfuzji i produkcji 
na parujących kropelkach. Równanie transportu uwzględniające te trzy czynniki wygląda 


następująco: 


—le;TJ+V(p,YV)=-V-J,+S,, (6.10) 


S, jest tutaj źródłem pary wodnej (parująca kropelka). 
Równanie ruchu kropelki 
Trajektorie kropelek były liczone przez całkowanie bilansu sił działających na kropelkę 


w Lagrange'owskim układzie odniesienia. W obliczeniach numerycznych uproszczono równa- 
nie 1.12 przez ograniczenie do wpływu siły Stokesa i grawitacji/wyporu. 
du 
mą = ŚTUR(u-v)+(pi=p,)V 8. (6.11) 
gdzie u jest prędkością kropelki, v - prędkością powietrza w pobliżu kropelki, u - lepkością 
dynamiczną powietrza, R — promieniem kropelki, P, - gęstością kropelki, P; - gęstością 


powietrza, V4 - objętością kropelki, g — wektorem przyspieszenia ziemskiego. 
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Parowanie 


Rozważany model uwzględnia wymianę masy przez parowanie wody z powierzchni 
kropelki. Podczas parowania, strumień pary na powierzchni kropelki zależy od gradientu 
koncentracji tej pary (model Stefana), w przybliżeniu jest to różnica pomiędzy koncentracją 


pary w otoczeniu kropelki i koncentracją pary przy powierzchni kropelki: 

NOSE AC, SC ah: (6.12) 
gdzie N, (kg mol/sm’) jest molekularnym strumieniem pary wodnej, k. jest współczynni- 
kiem transferu masy (m/s), Ć,, (kg mol/m’) — koncentracja pary na powierzchni kropelki, 


C,.„ (kg mol/m”) — koncentracja pary w otoczeniu kropelki. 


Koncentracja pary na powierzchni może być znaleziona z równania stanu gazu doskonałe- 
go przy założeniu, że przy powierzchni para jest nasycona a temperatura fazy gazowej jest 


równa temperaturze kropelki: 


DAT 4) 
= NO OE: 
C., TA (6.13) 


gdzie P,a jest ciśnieniem pary nasyconej (jest to funkcja temperatury), 7, jest temperatura 
kropelki. 
Ciśnienie pary nasyconej liczone jest z równania Clausiusa-Clapeyrona: 


L 


R, 


1 1 


Psal T)= Psat, oD Ti T 


(6.14) 


gdzie Pya.o= Psal To) a R, jest stałą gazową pary wodnej. 
Współczynnik transferu masy wyraża się liczbą Sherwooda i może być wyznaczony 
z równania [41]: 
= 2Rk c 1/20 1/3 
Sme a Sc `, (6.15) 
gdzie Sc jest liczbą Schmidta - Sc=u/(p,D,),a D, współczynnikiem dyfuzji pary wodnej. 
Zmiana masy kropelki jest liczona w modelu z równania: 


dm, 


dt 


=-N, HA, (6.16) 


gdzie H, jest masą molową pary wodnej Æa jest polem powierzchni kropelki. 
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Zmiana temperatury kropelki 
Kropelka wody unosząca się w powietrzu zmienia swą temperaturę wskutek wymiany 
ciepła z otoczeniem i parowania: 


dm, 


AL, (6.17) 


gdzie A4 jest powierzchnią całkowitą kropelki, h - współczynnikiem wymiany ciepła, L — 
ciepłem utajonym (ciepłem przemiany fazowej). Współczynnik h jest oszacowywany z kore- 


lacji Ranza — Marshalla [41]: 


h2R 1/2 1/3 
N,= =20+0.6Re, Pr, (6.18) 


u 


gdzie k jest współczynnikiem przewodnictwa cieplnego dla powietrza, Re, — liczbą Reynoldsa 


l 2R,|u— . ; 
dla kropelki ( Re,= paR ), Pr - liczbą Prandtla dla powietrza ( €„,% ). 

W symulacjach wykorzystany został program FLUENT 6.3.26 (Ansys Inc) [42]. Oprogra- 
mowanie to pozwala na rozwiązywanie metodą objętości skończonych przepływów z fazą 
dyskretną, czyli przepływów, w których występują cząsteczki ciekłe bądź stałe oddziaływają- 
ce z płynnym medium. 

Do dyspozycji są następujące opcje: 

e obliczenia trajektorii cząsteczek używając sformułowania Lagrange'a co pozwala na 

uwzględnienie bezwładności, oporu hydrodynamicznego, grawitacji; 
e parowanie cząsteczek ciekłych i wymianę masy z otoczeniem; 

e dyfuzję parujących składników; 

e wpływ energii przemian fazowych. 

Przy modelowaniu powietrza użyto opcji pozwalającej na uwzględnienie transportu skład- 
ników. W tym przypadku była to para wodna, której stężenie mogło fluktuować w czasie 
i przestrzeni. Modelowane powietrze składało się z dwóch substancji: suchego powietrza 
i pary wodnej. Ponieważ w każdym przypadku obliczeniowym ustalana była inna wilgotność 
początkowa (30%, 60%, 100%), różne było również stężenie pary wodnej. Dodatkowym 
źródłem pary wodnej były kropelki wody, które parowały w trakcie eksperymentu numerycz- 


nego. 


96 6 Badania numeryczne 


Obliczenia wykonano dla siatki obliczeniowej złożonej z 23550 elementów i geometrii 
płaskiej. W eksperymencie badane było oddziaływanie kropelek wody z zadaną strukturą 
wirową. Zasymulowany został pojedynczy dwuwymiarowy wir o rozkładzie prędkości opisa- 


nym następującymi równaniami: 


1 r” 
2 20° 


gdzie V,., V, oznaczają odpowiednio składową radialną i transwersalną pola prędkości 


(6.19) 


I 


V;=—r g AP 


w biegunowym układzie współrzędnych, r — odległość od środka wiru, o — odległość maksi- 


mum modułu prędkości od środka wiru, A - maksimum modułu prędkości (patrz rys.6.1). 


W przeprowadzonych symulacjach, maksimum prędkości wiru znajdowało się w odległo- 
ści lmm od środku wiru i wynosiło 2cm/s. Parametry te dobrano tak, aby były nieco większe 
od skal Kołmogorowa wynikających z danych doświadczalnych, aby w ten sposób sprawdzić 
jak wyglądają oddziaływania parujących kropelek z małymi skalami pola prędkości. Domeną 
obliczeniową było koło o promieniu 5mm. Ścianki ograniczające domenę były adiabatyczne 


z kinematycznym warunkiem poślizgu na brzegu domeny obliczeniowej. 


r [mm] 


Rys. 6.1: Profil prędkości transwersalnej na przekroju przez wir dla warunku początkowego. 
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y [mm] 


Rys. 6.2: Warunki początkowe dla modelu wiru: izolinie modułu prędkości i początkowe położenie 12 
kropelek. Pole grawitacji skierowane przeciwnie do zwrotu osi Y. 

Wyniki symulacji przedstawione są na wykresach (rys. 6.3). Symulacje dla wilgotności 
100% potraktowane zostały jako referencyjne (w tych warunkach nie ma parowania). Przed- 
stawione wyniki pokazują różnice energii kinetycznej pomiędzy symulacją referencyjną 
a symulacjami dla wilgotności 30% i 60%. W obu przypadkach widać, że zmianie ulega 
głównie pionowa składowa pola prędkości (por. rys. 6.2 i rys. 6.3). Z rys 6.2 wynika, że po 
lewej stronie od centrum wiru prędkość jest skierowana przeciwnie do grawitacji, natomiast 
po prawej stronie zgodnie z grawitacją. Można zauważyć, że energia kinetyczna jest większa 
w obszarach, gdzie pole prędkości jest dodatnio skorelowane z grawitacją, natomiast energia 


kinetyczna jest mniejsza w obszarach, gdzie ta korelacja jest ujemna. 


Na podstawie tych symulacji można stwierdzić, że parowanie kropelek oddziałuje bezpo- 
średnio tylko na pionowe fluktuacje pola prędkości. Pionowe ruchy opadające zostają wzmoc- 
nione wskutek oddziaływania hydrodynamicznego kropelek i zwiększonej gęstości powietrza 
spowodowanej parowaniem. Pionowe ruchy wstępujące poddane tym samym mechanizmom 
ulegają wygaszeniu. Wielkość tego efektu wskazuje, że proces parowania jest czynnikiem 


odpowiedzialnym za istotny wzrost energii fluktuacji pola przepływu. 


Połączenie przedstawionego modelu z symulacją przepływu w całej komorze chmurowej 
aczkolwiek możliwe, nie zostało zrealizowane z uwagi na ograniczone moce obliczeniowe 
w Zakładzie. Prace takie są m. in. prowadzone we współpracy z S. Malinowskim, M. Andrej- 


czukiem, W.W. Grabowskim, P.K. Smolarkiewiczem [2] 
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t=0.02s 


t=0.04s 


t=0.06s 


t=0.08s 


t=0.10s 


x [mm] 


Rys. 6.3: Porównanie symulacji numerycznych dla wilgotności 30% - lewa kolumna i wilgotności 60% 
- prawa kolumna. Kontury prezentują różnicę energii kinetycznej w stosunku do symulacji dla 
wilgotności 100% (brak parowania). Pole grawitacji skierowane przeciwnie do zwrotu osi Y. 


7 Podsumowanie 


Przedmiotem przeprowadzonych badań była drobnoskalowa turbulencja w chmurach. 
Wiedza w tym temacie jest nadal niewystarczająca do wyjaśnienia pewnych zjawisk obserwo- 
wanych w atmosferze. Głównym powodem jest brak danych eksperymentalnych o rozkładach 
prędkości kropelek chmurowych. Celem pracy było zapełnienie tej luki przez zastosowanie 
w tej dziedzinie anemometrii obrazowej i laboratoryjnego modelu chmury. Zbadano proces 
mieszania się chmury z czystym powietrzem — próbując odtworzyć zjawiska zachodzące na 
granicy chmury kroplowej i powietrza w rzeczywistej chmurze. Przyjęto w pracy założenie, 
że główną rolę w procesach generacji energii turbulencji odgrywa proces mieszania na grani- 
cy chmury kroplowej i czystego powietrza. Takie założenie uzasadnia stopień filamentacji 
struktur chmurowych obserwowany w warunkach atmosferycznych i w komorze laboratoryj- 
nej, wskutek którego efektywna powierzchnia mieszania przekracza o kilka rzędów wielkości 


powierzchnie objętości zajmowanej przez strukturę chmurową. 


Zastosowana w pracy metoda anemometrii obrazowej (PIV) pozwoliła na zbadanie turbu- 
lentnych charakterystyk pól prędkości z rozdzielczością poniżej skali Kołmogorowa. Oznacza 
to, że odtworzono strukturę przepływu z dokładnością konieczną dla pełnego opisu przepływu 
turbulentnego, umożliwiając tym samym porównania z symulacjami numerycznymi realizo- 
wanymi metodami DNS dla małych skal turbulencji. W ramach badań eksperymentalnych 
wyznaczono między innymi zależność współczynnika dyssypacji lepkiej w funkcji wilgotno- 
ści względnej i od stopnia zmieszania (odległości od wlotu do głównej komory). Stwierdzono 
jakościową zgodność tych badań z zaproponowanym modelem mieszania uwzględniającym 
zmiany gęstości mieszaniny w funkcji wilgotności i proporcji zmieszania, związane z prze- 
mianami termodynamicznymi. Potwierdza to hipotezę o wpływie przemian fazowych na 


strukturę turbulencji w małych skalach przepływów atmosferycznych. 


W badaniach turbulentnych charakterystyk pól przepływu wykorzystano funkcje struktury 
wyznaczone dla danych doświadczalnych. Stwierdzono istotny wpływ wilgotności i proporcji 
zmieszania na badaną tą metodą strukturę turbulencji. Dodatkowym argumentem wskazują- 
cym na wpływ parowania kropel na strukturę lokalnej turbulencji w chmurze są eksperymenty 
przeprowadzone dla zawiesiny kroplowej złożonej z nieparującej cieczy. Wyraźna zmiana 


funkcji struktury i współczynnika dyssypacji lepkiej potwierdza, że parowanie kropel jest 
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istotnym czynnikiem wpływającym na przepływ. Przeprowadzone dodatkowo badania nume- 
ryczne parowania pojedynczej kropli potwierdziły, że obecność parujących kropelek w struk- 


turze wirowej ma bezpośredni wpływ na ewolucję pola energii kinetycznej. 


Wszystkie wymienione powyżej wyniki wskazują, że parowanie w wyniku niejednorodne- 
go mieszania mas powietrza ma istotny wpływ na turbulencję. O ile wpływ przemian fazo- 
wych na wielkoskalowe ruchy konwekcyjne w atmosferze, które są źródłem turbulencji jest 
od dawna oczywisty, o tyle wpływ parowania w małych skalach nie jest nadal do końca 
poznany. Na podstawie wykonanej pracy można przypuszczać, że wpływ ten jest istotny. Jest 
to tym bardziej ważne, że właśnie w małych skalach zachodzą procesy koalescencji kropelek 
co prowadzi do zmiany widma wielkości kropelek i może prowadzić do takich procesów jak 


ciepły deszcz. 


Funkcje struktury jak i statystyki pola prędkości pokazują anizotropię turbulencji. Jak 
pokazano na podstawie funkcji struktury, anizotropia ta dotyczy szerokiego zakresu skal łącz- 
nie ze skalą Kołmogorowa. Statystycznie fluktuacje pionowe są silniejsze od fluktuacji pozio- 


mych. 


Przeprowadzone badania dostarczyły cennych danych eksperymentalnych do porównań 
z modelami numerycznymi. Wstępne symulacje numeryczne wykonane techniką DNS dla 
małych skal turbulencji [2] potwierdziły obserwowany w eksperymencie wpływ parowania 
i fluktuacji wypornościowych na strukturę pola fluktuacji turbulentnych. Kolejnym koniecz- 
nym krokiem jest obecnie powiązanie tych efektów ze statystyką zderzeń kropel, a więc 
znalezienie jaki wpływ mają te efekty na intensyfikację koalescencji kropel w rzeczywistej 


chmurze. 
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Mosix Cluster [43]. 


9 Dodatki 


9.1 Obliczenia diagramu mieszania 


Rozważmy proces mieszania nasyconego powietrza zawierającego wodę kropelkową 
z powietrzem nie zawierającym wody kropelkowej nazwanym tutaj powietrzem czystym. 
Mieszanie to odbywa się bez wymiany masy i energii z otoczeniem i przy zachowaniu stałego 
ciśnienia. Powietrzem suchym będziemy nazywać powietrze bez pary wodnej. Dowolną masę 
powietrza można podzielić na powietrze suche i parę modną. Wprowadźmy następujące ozna- 


czenia: 
mai — masa powietrza suchego zawartego w pierwszej cząstce, 
Ma2 — masa powietrza suchego zawartego w drugiej czastce, 


Ma — masa powietrza suchego po mieszaniu, 


mM, — masa pary wodnej zawartej w pierwszej cząstce, 
m,» — masa pary wodnej zawartej w drugiej cząstce, 


m, - masa pary wodnej po zmieszaniu, 


T, — temperatura pierwszej cząstki, 
T, — temperatura drugiej cząstki, 
Tn — temperatura po zmieszaniu bez parowania , 


T — temperatura końcowa, 


X= - - stosunek zmieszania pary wodnej dla pierwszej porcji powietrza, 


X= mę - stosunek zmieszania pary wodnej dla drugiej porcji powietrza, 
a2 
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m 
X =— - końcowy stosunek zmieszania pary wodnej, 
m 


a 


m 


wl > : $ $ > 
w= mda stosunek zmieszania wody kropelkowej w pierwszej cząstce, 
al 
= M . 5 z EZ, 
kz stosunek zmieszania wody kropelkowej w drugiej cząstce, 
a2 


m 


w 


w= = końcowy stosunek zmieszania wody kropelkowej, 


R,=461 E - stała gazowa pary wodnej, 


a 


z =0.622 -stosunek stałej gazowej czystego powietrza do stałej gazowej pary wodnej. 


v 


Zakładamy, że znamy parametry opisujące obie cząstki, takie jak: temperatura (7), 72), 
stosunek zmieszania pary wodnej (X; i X2), stosunek zmieszania wody kropelkowej (w; i w2) 
i ciśnienie, które jest równe ciśnieniu atmosferycznemu. Zastanówmy się jaka będzie tempe- 
ratura (T), stosunek zmieszania pary wodnej (X) i stosunek zmieszania wody kropelkowej (w) 


po zmieszaniu obu porcji powietrza w zadanych proporcjach zmieszania: 


M, M, ma ma 
k —_—* =, k= —=— (9.1) 
Mi EM Ma Mt Ma2 Ma 
Rozbijmy rozważany proces na dwa etapy: 
e W pierwszym następuje tylko mieszanie (bez parowania wody). Obliczymy 
najpierw jaka będzie temperatura (7T,), i stosunki zmieszania (X, 1 Wn) takiej mieszan- 
ki. 
e W drugim etapie obliczymy zmianę tych wielkości (AT, AX, Aw) po odparowaniu 
wody kropelkowej. Uwzględnimy przy tym fakt, że parowanie będzie następować do 
chwili, kiedy ciśnienie pary wodnej osiągnie stan nasycenia lub gdy wcześniej wypa- 


ruje cała woda kropelkowa zawarta w rozważanym przez nas układzie. 


Gdy nie ma parowania to z bilansu energii otrzymujemy następujące równanie (na tym 
etapie pomijamy tutaj zmiany objętości): 


(m, C pat My C py tM, Gl he. tm cto.) Ta (9 2) 
=( M,C atM, C „A Myć „JB, j 


w pw n 
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Ponieważ masa wody i pary wodnej jest mała w stosunku do masy powietrza suchego może- 
my je zaniedbać, dodatkowo dzieląc równanie 9.1 przez c,,m, otrzymujemy: 

k Titk, =f. (9.3) 
W pierwszym etapie nie zmienia się całkowita masa wody kropelkowej i masa pary wodnej. 
Z bilansu masy wody kropelkowej łatwo więc pokazać, że stosunek zmieszania ciekłej wody 
będzie wynosić: 

w,„=k,w;,tk,w,. (9.4) 
W naszym przypadku tylko jedna porcja powietrza zawiera wodę kropelkową, zatem: 

w„=k,w,. (9.5) 

Podobnie z bilansu masy pary wodnej wynika, ze: 

X „~k X tk, X. (9.6) 
Znamy już więc parametry termodynamiczne naszego układu po zmieszaniu bez parowania 


wody. Teraz włączamy parowanie. 


Parowanie spowoduje zmianę temperatury, zmianę zawartości wody kropelkowej i pary 
wodnej. Zmianę temperatury można policzyć z bilansu energii: 

AU=Q+W. (9.7) 
Praca w tym przypadku jest związana ze zmianą objętości w stałym ciśnieniu. Przyjmując, że 
wpływ zmiany objętości jest zaniedbywalnie mały to zmiana energii wewnętrznej jest spowo- 


dowana pobraniem z układu ciepła parowania 


m,C,,AT=LAm,,, (9.8) 
skąd otrzymujemy dalej: 
L 
AT=—A 
Bao (9.9) 


pa 
Temperature jaka ustali się po parowaniu można zapisać równaniem: 


T=T„+AT=k,T;+k, T+ Aw, (9.10) 
Aby obliczyć te temperature musimy znać jaka i wody wyparuje. Woda paruje do osią- 
gnięcia stanu nasycenia. Wstawmy więc w równanie 9.8 zależność od ciśnienia pary nasyco- 
nej. Na początku zamieńmy stosunek zmieszania ciekłej wody na stosunek zmieszania pary 
wodnej: 

Po o Pa (9.11) 
Cha 

co wynika z zachowania całkowitej masy wody w układzie. 
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AX=X(T)-X,„. (9.12) 
Temperatura T jest w tym przypadku temperaturą przy której powietrze jest nasycone, zatem 
X(T) jest stosunkiem zmieszania pary wodnej dla temperatury nasycenia. Wykorzystując 
zależność między ciśnieniem pary wodnej a stosunkiem zmieszania pary wodnej i wstawiając 
do niej prężność pary nasyconej E,(T)’ otrzymujemy: 
EE (T 
a 
p-E;(T) 
gdzie p jest ciśnieniem, a €=R,/R,=0.622 | 


(9.13) 


Ponieważ ciśnienie pary wodnej jest małe w stosunku do ciśnienia atmosferycznego może- 


my uprościć powyższą zależność pozbywając się ciśnienia pary wodnej z mianownika: 


X(T)="ZEs(T), (9.14) 
gdzie ciśnienie pary nasyconej otrzymujemy z równania Clausiusa-Clapeyrona: 
L|/1 1 
E;(T)=E |= 
s(T)=E exp RAT, T|) (9.15) 


Gdzie E;„=E;(T,) jest znaną wartością ciśnienia pary nasyconej zmierzonej w temperatu- 


rze To. 
Z równań 9.2, 9.4, 9.9, 9.10 otrzymujemy: 


L 
T=k Tikt |E EST) X-k). (9.16) 


pa 


Korzystając dalej z faktu, że k,+k,=l możemy z powyższej równości pozbyć się K,: 


L 
T=(1-k;) Pitka Z$ E(T)-0-4) Xba 


pa 


; (9.17) 


Aby znaleźć szukaną zależność T(k,) należy powyższe równanie rozwiązać numerycz- 


nie. 


Wyprowadzone powyżej równanie opisuje układ, w którym osiągnięty jest stan nasycenia 
parą wodną. W badanym procesie sytuacja taka będzie obserwowana tylko wtedy, gdy 
w układzie jest wystarczająco dużo wody kropelkowej. Jeśli wody kropelkowej jest mniej to 


wyparowuje ona cała. 
Mamy więc dwa reżimy: 


e dla k „E|0. k „| wyparowuje cała woda — ten przypadek opisuje równanie 


I Ciśnienie pary nasyconej £.(T) zależy od temperatury i wiąże się z definicją wilgotności względnej. Jeśli 
e(T) jest aktualnym ciśnieniem pary wodnej, to wilgotność względna wyraża stosunek aktualnego ciśnienia 
do ciśnienia pary nasyconej w danej temperaturze: RH =100%:e ( T )/ E; ( T ) ; 
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L 
T=(1-k,)Ti+k,T, =-= Wyk, (9.18) 
pa 


e dla k Elk, 1| osiągnięty jest stan nasycenia po odparowaniu odpowiedniej masy 
wody (część wody kropelkowej pozostaje w postaci ciekłej) — ten przypadek opisuje 


równanie 9.17. 


kK, jest tutaj wartością graniczną dla której po odparowaniu całej wody powietrze jest nasy- 


cone. 


Zależność T(k) pozwala nam wyznaczyć zmiany temperatury gęstościowej T, w funkcji 


stosunku zmieszania przedstawione na rys. 2.1. 


9.2 Opis algorytmu poszukiwania lokalnych przemieszczeń 


Cyfrowe obrazy można traktować jak tabele, czy macierze. Indeksy poszczególnych 
elementów macierzy są związane liniowym przekształceniem ze współrzędnymi na płasz- 
czyźnie fotografii i będą tutaj utożsamiane jako współrzędne - odpowiednio indeks i — współ- 
rzędna pozioma, j - pionowa. Wartości poszczególnych elementów macierzy odpowiadają 
poziomowi jasności w danym punkcie fotografii. Jasność przybiera całkowite wartości dodat- 
nie z zakresu określonego przez właściwości kamery i liczbę bitów przetwornika analogowo- 


-cyfrowego. 


Niech $7, $> oznaczają odpowiednio macierze odpowiadające cyfrowym obrazom o tych 
samych rozmiarach MxN. S;(ij), Sij) niech oznaczają odpowiednio elementy tych obrazów 
(piksele) o współrzędnych i,j. Współczynnik korelacji można zatem zdefiniować w następują- 


cy sposób: 


M= 
Mz 


(Sili, f)—(Sy)) (S214, J)—45,)) 


Cc(8,8,)=——— — , (9.19) 
|> 2 (5,6080) 2,2 SNES 
i=1 j=l i=l j=l 
gdzie symbol <> oznacza uśrednianie według następującego wzoru: 
1 N 
Sp =—— S (i,j). 52 
( » wan > (i j) (9.20) 


Można pokazać, że w liczniku wyrażenia z równania 9.19 prawdziwa jest następująca 


zależność: 
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(Sili -SSi 1)-(S))=2, 2 (S,(i,5)-(80)(5(i,J)-Z), (9.21) 


i=1 j=l 


Mes 
Maz 


Il 
jeż 

II 
jak 


i=1j 
gdzie Z jest dowolną stałą. Wstawienie Z=0 pozwala na uproszczenie i przyspieszenie obli- 


czeń. Podobnie można uprościć występujące w mianowniku 9.19 dyspersje: 


Ms 


(Si, J-(S)P= DD (SiG, JP -MN (Sy). (9.22) 


i=1 j=l 


II 
= 

II 
— 


i=1j 


Uproszczenia te są szczególne przydatne, jeśli porównujemy jeden obraz ze zbiorem obra- 


ZOW. 
Załóżmy, że mamy ustaloną próbkę S; i chcemy porównać ją z wieloma różnymi próbkami S>. 


Dla próbki $, liczymy tylko raz macierz Š, według wzoru: 
Ś,(i,j)=S,(i, j)-(S,) (9.23) 


i dyspersję macierzy $7: 


o- Sisa DSDHA HEA. (9.24) 


Wyrażenie na współczynnik korelacji 9.19 możemy zatem przekształcić do następującego 


wzoru: 


CSS) >. (9.25) 
a 


Znajdowanie podobnych fragmentów na dwóch fotografiach opiera się na sprawdzaniu podo- 


bieństwa małego fragmentu z pierwszej fotografii ze zbiorem fragmentów z drugiej fotografii. 


Niech O i S będą obrazami cyfrowymi. Dla uproszczenia załóżmy ze są to obrazy kwadra- 
towe. Niech O ma zatem wymiar LxL a S wymiar MxM, przy czym wymiar L>M. Operację 
korelacji wzajemnych wykonuje się, aby znaleźć wewnątrz obrazu O fragment najbardziej 
podobny do obrazu $. W tym celu korzystając ze współczynnika korelacji porównuje się 
wszystkie fragmenty o wymiarach MxM zawarte wewnątrz obrazu O z obrazem S. Otrzymu- 
jemy zatem nie jedną wartość współczynnika korelacji, ale zależność tego współczynnika od 
położenia testowanych fragmentów. Zależność tę można opisać wzorem: 

M M 

2 2Śli.J)O(i+x,j+y) 

_ i=l j=l 4 (9.26) 
MAKOŻ)-(0,„)* 
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gdzie 
Oa O(x+i, y+ j) (9.27) 
(0 )= ar È È Olti y+) (9.28) 


Tak zapisana funkcja korelacji jest szczególnie wygodna dla zastosowania standardowych 
algorytmów przetwarzania obrazów z wykorzystaniem transformaty Fouriera. Iloczyn w licz- 
niku można traktować jako splot dwóch periodycznych funkcji, z których jedna — zawierająca 
obraz próbki ma ograniczony nośnik o rozmiarze próbki (to znaczy sygnał jest ograniczony 
tylko do kwadratu MxM, a reszta przestrzeni jest wypełniona zerami). Na mocy twierdzenia 
o splocie transformatę Fouriera splotu można wyrazić jako iloczyn transformat Fouriera spla- 


tanych funkcji: 
F|f0g)=F|F|F|g]. (9.29) 


gdzie F oznacza transformatę Fouriera, a symbol © oznacza operację splotu. 
Odpowiada to wyrażeniu dla korelacji: 

F|f*g]=F|f)F(s). (9.30) 
gdzie symbol % oznacza operację korelacji wzajemnej (z ang. cross correlation), a symbol * 


- sprzężenie. Zastosowanie do tych obliczeń szybkiej transformaty Fouriera (FFT) pozwala na 


znaczne przyspieszenie działania programu. 


Ponadto wyrazy występujące w mianowniku <O0,„> i <O’,,> nie zależą od próbki S. Nie 
trzeba zatem liczyć ich za każdym razem, gdy wybierana jest kolejna próbka S. W programie 
wyrazy te liczone są na początku dla całej fotografii jeszcze przed przystąpieniem do liczenia 


lokalnych korelacji wzajemnych. 


W przypadku opracowanego przez autora programu PIV-Kor próbką jest kwadratowy frag- 
ment pierwszej fotografii natomiast obszarem poszukiwań kwadratowy obszar drugiej foto- 
grafii w którym chcemy znaleźć fragment mu odpowiadający. Wielkość tego obszaru jest 
uzależniona od maksymalnych przemieszczeń, których się spodziewamy w przepływie. Przy 
czym współrzędne środka obu kwadratów są takie same. Wtedy C(x,y) można interpretować 
jako zależność współczynnika korelacji pomiędzy ustalonym fragmentem z pierwszej fotogra- 
fii, a fragmentem z drugiej fotografii, który odpowiada hipotetycznemu przemieszczeniu 


próbki o wektor (x,y) z początkowego położenia. C(x,y) wyznacza prawdopodobieństwo prze- 
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mieszczenia (x,y). Znalezienie maksimum tego prawdopodobieństwa pozwala na oszacowanie 


przemieszczenia fragmentu S$ i prędkości grupy cząsteczek zawartej w tym fragmencie. 


9.3 Dyskusja wiarygodności techniki PIV 


Opisany w niniejszej pracy program PIV-Kor stosowany do pomiarów PIV w niniejszych 
badaniach opiera się na metodach statystycznych, które pozwalają na oszacowanie prawdopo- 
dobieństwa przemieszczeń fragmentów obrazów. Dla poprawienia dokładności wyznaczania 
przemieszczeń zastosowano w nim również kilka dodatkowych usprawnień jak deformacje 
obrazów i dynamiczne zmniejszanie fragmentów w kolejnych iteracjach, filtrowanie. Zacho- 
dzi więc pytanie na ile zastosowane algorytmy pozwalają na zwiększenie dokładności pomia- 
ru, a na ile same przyczyniają się do powstania artefaktów? Na ile można uzyskane wyniki 
obliczeń traktować jako przemieszczenia fragmentów obrazu fizycznego przepływu? Może 
zachodzić podejrzenie, że różne algorytmy PIV mogą dawać zupełnie inne wyniki. W tej 
części pracy porównano wyniki działania napisanego na potrzeby niniejszej pracy programu 
PIV-Kor i komercyjnego oprogramowania VidPIV firmy ILA GmBh. Do porównania użyto 
pary fotografii z komory chmurowej. Wybrana fotografia charakteryzuje się dużą niejedno- 


rodnością w rozkładzie przestrzennym kropelek (patrz rys. 9.1) 


a 
Q 
© 
T 
i 


y [piksele] 


‘ 
a 


900 |- Poe E J 


100 300 500 700 300 1100 
x [piksele] 

Rys. 9.1: Pierwsza fotografia z pary użytej do porównania oprogramowania PIV. Ciemne obszary 

odpowiadają dużej koncentracji kropelek. Jasne obszary odpowiadają obszarom o małej koncentracji. 

Nawet obszary widoczne jako zupełnie białe zawierają obrazy kropelek, które są wykrywane przez oba 

programy, dzięki czemu pole prędkości mogło być oszacowane na całej powierzchni fotografii. 
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Rys. 9.2: Porównanie pól przemieszczeń otrzymanych z obu programów. Pierwsza kolumna — 
składowa pozioma przemieszczenia, kolumna druga — składowa pionowa przemieszczenia. Górny rząd 
— wyniki dla programu VidPIV firmy ILA, dolny rząd — wyniki dla własnego oprogramowania — PIV- 
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Rys. 9.3: Porównanie profili przemieszczenia wzdłuż poziomych przekrojów przez pole prędkości. 
Lewa kolumna — składowa pozioma przemieszczenia, prawa kolumna — składowa pionowa 


przemieszczenia. Linia ciągła — 


VidPIV, linia przerywana — PIV-Kor. 
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Rys. 9.4: Porównanie profili przemieszczenia wzdłuż pionowych przekrojów przez pole prędkości. 
Lewa kolumna — składowa pozioma przemieszczenia, prawa kolumna — składowa pionowa 
przemieszczenia. Linia ciągła — VidPIV, linia przerywana — PIV-Kor. 


Z analizy rys 9.2-9.4 wynika, że obydwa programy generują zgodne wyniki dla prawie 
całej powierzchni fotografii. Różnica pojawia się w niewielkim fragmencie u dołu fotografii 
(x~1000 pikseli, y~900 pikseli). Nie udało się ustalić, który z programów generuje w tym 
fragmencie błędne wyniki. Problem polega raczej na tym, że w tym obszarze pojawiają się 
prawdopodobnie duże fluktuacje składowej prędkości w kierunku prostopadłym do fotografii, 
co powoduje w skrajnych przypadkach utratę informacji o przemieszczeniach. Odpowiednie 
dobranie interwału czasowego pomiędzy obydwiema fotografiami pozwala na zmarginalizo- 
wanie tego efektu, ale w trójwymiarowych, silnie turbulentnych przepływach całkowite 
wyeliminowanie go jest raczej niemożliwe. Problem ten pojawia się jednak na niewielu foto- 


grafiach i dotyczy stosunkowo niewielkiej powierzchni objętej pomiarem. 


Tak dobra zgodność w działaniu dwóch różnych i niezależnie skonstruowanych progra- 
mów pozwala przypuszczać, że pomiary przemieszczeń fragmentów obrazów wykonane tech- 


niką PIV dla zdjęć z komory chmurowej są wiarygodne. 


Jak to zostało opisane wcześniej, program PIV-Kor nie mierzy przemieszczeń pojedyn- 
czych cząsteczek znacznika, a przemieszczenie grupy cząsteczek znajdujących się na niewiel- 
kim fragmencie obrazu. Pomiar ruchu pojedynczych kropelek nie był możliwy, ze względu na 
bardzo nierównomierny rozkład przestrzenny kropelek i występowanie obszarów, w których 


było ich za dużo na to, żeby odseparować obrazy pojedynczych kropelek. Zmierzona pręd- 
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kość w danym punkcie jest więc uśrednioną prędkością grupy kropelek znajdujących się 
w otoczeniu tego punktu. Aby zmniejszyć efekt tego uśredniania, otoczenie punktu, czyli 
wielkość fragmentu próbki obrazu, dla którego szukane jest przemieszczenie jest zmniejszana 
w kolejnych iteracjach działania programu. W połączeniu z deformacją obrazu pozwala na 
dokładniejsze wyznaczenie przemieszczeń, jednak pozostaje pytanie na ile taki algorytm ma 
wpływ na pomiar. 

Aby zbadać wpływ krzywizny trajektorii i efektu uśredniania prędkości, program PIV-Kor 
został przetestowany na sztucznej parze fotografii testowych wygenerowanych numerycznie 
dla zadanego, znanego pola prędkości [44]. Dzięki temu możliwe jest po zastosowaniu 
programu PIV do takiej pary obrazów, porównanie zmierzonego pola prędkości z chwilowym 
polem prędkości znanego przepływu użytego do wygenerowania tych fotografii. Porównanie 
takie pozwoli sprawdzić, jak niedoskonałości zastosowanej techniki PIV wpływają na dokład- 


ność pomiaru. 


W pomiarach użyto przepływu wokół dwóch walców (rys. 9.5 i 9.6). Średnie przesunięcie 
dla oryginalnego pola wynosiło w tym przypadku 7,6 piksela. Obliczono różnice pomiędzy 
wektorami pola oryginalnego a wektorami uzyskanymi z pomiarów. Średni moduł tej różnicy 
uzyskany przez uśrednianie po całej powierzchni pomiaru wynosił około 0,3 piksela przy 
maksymalnej różnicy około 0,8 piksela. Wynika z tego, że pomiary PIV z wykorzystaniem 
programu PIV-Kor odtwarzają pole prędkości z dokładnością rzędu 0,3 piksela co dla typo- 


wych przemieszczeń oznacza błąd względny wektora prędkości poniżej 3%. 


Powyższa dyskusja pokazuje, że pomimo znanych niedoskonałości, odpowiednio stosowa- 
na technika PIV wraz z napisanym na potrzeby niniejszej pracy oprogramowaniem PIV-Kor, 
umożliwia uzyskanie wiarygodnych pomiarów pola prędkości kropelek chmurowych 


w warunkach prezentowanego tutaj eksperymentu. 


Rys. 9.5:Pierwszy z dwóch obrazów użytych do testowania programu PIV-Kor. 
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Rys. 9.7: Pole błędu w pikselach — moduł różnicy między wektorami pola oryginalnego a wektorami 
pola wyznaczonego programem PIV-Kor. Widoczne jest, że najwięcej problemów algorytm miał w 
obszarach dużego ścinania prędkości przy powierzchni walców. 


9.4 Tabela z wynikami pomiarów wewnątrz komory 


Poniżej prezentowana jest tabela wyników z 50 serii pomiarowych przeprowadzonych 
wewnątrz komory chmurowej z udziałem kropelek wody. Podczas badań wykonano również 
wiele eksperymentów próbnych i pomocniczych, których wyniki nie są w tej pracy publiko- 


wane. 
Zastosowano następujące oznaczenia: 
- Nr- numer eksperymentu w kolejności chronologicznej, 
- Poz. - odległość pomiaru od wlotu do głównej komory, 
- RH- wilgotność względna czystego powietrza wewnątrz komory, 
-  €— współczynnik dyssypacji lepkiej (5.7), 


- n- skala długości Kotmogorowa (5.10), 
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(u), (w) - uśredniona składowa pozioma i pionowa prędkości, 


0,, O, - odchylenie standardowe dla składowej poziomej i pionowej (5.11), 


S(u), S(w) — skośność dla składowej poziomej i pionowej pola prędkości (5.12), 
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K(u), K(w) — spłaszczenie dla składowej poziomej i pionowej pola prędkości (5.13), 


AM, , A„ - mikroskala Taylora dla składowej poziomej i pionowej pola prędkości (5.2). 


Tabela 8: Zestawienie wyników pomiarów właściwości przepływu wewnątrz komory z udzia- 
łem kropelek wody. Odpowiadające im pola prędkości zostały zarchiwizowane na nośnikach 


DVD. 
RH e€ u w o T, ATES 

Nr |Poz. | poż] | [m/s | [min] a: MA Ewa EAC ENEA 
1| 70| 20| 0,016| 0,68 |-0,002 0,30| 0,05 0,062 | -0,04 | 0,03] 0,0 -0,2| 54) 6,6 

2) 70| 42| 0,007) 0,84 0,009 0,26 | 0,042 0,052 | -0,09 | 0,09 0,0 -02| 74| 93 

3) 70) 49| 0,006| 0,88 | 0,009) 0,24 | 0,044 0,052 | -0,08 | 0,01 -0,1 -0,2| 89 11,0 

4| 70| 61| 0,006) 0,89 | 0,013 0,23 | 0,038 0,052 | -0,03 0,01 0,0 -01| 80) 1,2 

5| 601 25| 0,011| 0,75 -0,015 0,28 0,048 0,059 | -0,02 | 0,07 |-0,1 -03| 62| 7,8 

6 60) 20| 0,012| 0,73) 0,001) 0,30 | 0,049 0,056  -0,10 | 0,01 -0,1 -01| 60) 7,0 

7| 60| 29| 0,010| 0,77 | -0,007 0,27 | 0,045 |0,056| 0,03 | 0,11 |-0,2 -04| 64| 7,8 

8 60| 40) 0,009| 0,79 |-0,001 | 0,28 | 0,046 0,054 | -0,06 |-0,05 | 0,3 -0,2| 6,8| 8,1 

9 60) 62| 0,005| 0,90) 0,004) 0,21 | 0,039 0,047 | -0,06 |-0,02 -0,1 -0,2| 8,8 11,0 

10 50, 62| 0,007| 0,84) 0,008 0,24 | 0,042 0,055 | -0,08 | 0,00 -0,1 -03| 7,5 10,0 
11 50 25 0,014 0,71| 0,013) 0,22 | 0,042 0,057 0,01 | 0,00 |-0,1 -0,1 70) 6,7 
12 50| 32| 0,012| 0,73) 0,002 0,27 | 0,046 |0,064 | -0,07 | 0,06 | 0,2 | -02| 5,8| 7,8 
13) 50| 42| 0,010| 0,77) 0,001 | 0,25 | 0,047 0,062! 0,01 |-0,04 |-0,2 -04| 65| 8,9 
14) 50| 50| 0,007| 0,85 |-0,001| 0,21 | 0,038 0,054 | -0,11| 0,17 0,1| 00| 72) 9,7 
15) 50| 61| 0,008| 0,81 | 0,003 0,22 0,045 0,056 | -0,12 |-0,05 -0,1 | -0,3 7,2| 93 
16) 40| 59| 0,009| 0,79) 0,014) 0,22 | 0,045 0,062 | -0,18 | 0,10 0,2 -0,5| 67) 9,2 
17) 40) 20| 0,014| 0,71| 0,002) 0,27| 0,05 0,064) 0,04 |-0,14 |-0,3 -03| 54| 69 
18) 40| 32| 0,012| 0,73) 0,000) 0,24 | 0,044 0,062 | -0,07 | 0,02 | 0,0 | -01| 56. 7,6 
19 40) 41| 0,011| 0,75) 0,002 0,25 0,046 0,062 | -0,07 | 0,00 | 0,0 -0,3 59| 7,8 
20| 40| 52| 0,010) 0,77| 0,016) 0,24 | 0,044 [0,063 | 0,06 | 0,06 -0,1 -03| 6,0] 9,0 
21| 401 46| 0,011) 0,75| 0,006 | 0,24 | 0,046 (0,057 | -0,08 | 0,03 -0,1 -03| 60) 7,4 
22 40| 48| 0,008! 0,80| 0,010) 0,21 | 0,04 0,056 | -0,09 | 0,15 -0,1 -02| 6,0] 8,2 
23| 30| 56| 0,011 | 0,75| 0,011| 0,16 | 0,043 |0,062 | -0,06 | 0,03 | 0,0 | -03| 5,5| 7,7 
24 30| 65| 0,007  0,83| 0,013) 0,18 | 0,037) 0,05 |-0,01 | 0,15| 0,0, 0,01 64| 83 
25| 30| 20| 0,021 | 0,64 |-0,001 | 0,21 | 0,053 | 0,064 | -0,07 |-0,13 |-0,1 |-0,1 | 46| 5,5 
26| 30| 36| 0,017| 0,67 | -0,011| 0,18 | 0,044 |0,057 |-0,08 | 0,04 | 0,3 |-0,2| 43| 5,5 
27| 30| 23| 0,021  0,64| 0,002) 0,19 | 0,052 0,067 | -0,02 |-0,16| 0,0 -01| 4,6] 6,0 
28  50| 23| 0,014  0,70| 0,006) 0,23 | 0,052 0,061 | 0,04 |-0,27| 0,1 -02| 43| 4,9 
29| 70| 23| 0,016) 0,69| 0,005) 0,27] 0,05) 0,06 | -0,04 | 0,00 |-0,1 | -0,2| 4,5| 5,4 
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RH © n lu) w) || o On Ape | 
[%] | [m’/s*] | [mm] | [m/s] | [m/s] | [m/s] | [m/s] | ©) | SY) ROD) ROW) 


w 


Nr | Poz. pera ee 


30 70 23) 0,015) 0,69) 0,008 | 0,25 | 0,048 | 0,058 | -0,10 | 0,13 | 0,2} 0,1) 5,2) 6,1 


31} 70, 23) 0,014) 0,71 | 0,013 | 0,25 | 0,046 | 0,052 | -0,01 | 0,02 | 0,1 | -0,1| 5,7) 6,3 


32} 70, 23| 0,014) 0,70 | -0,002 | 0,25 | 0,047 | 0,054 | -0,11 | 0,14 0,3) 0,2) 4,8) 5,6 


33} 70 23| 0,015) 0,69) 0,008 | 0,26 | 0,049 | 0,06 | -0,08 | 0,13) 0,1; 0,1) 5,2) 6,1 


34} 70, 23| 0,014) 0,71 | 0,000 | 0,26 | 0,046 | 0,053) 0,06 | 0,03 | 0,0| -0,2) 5,6) 6,4 


35) 70 38) 0,008) 0,82) 0,012) 0,24 | 0,043 | 0,059 | -0,11 | 0,07 | 0,0 | -0,2| 7,3) 10,0 


36} 50, 46) 0,007) 0,83) 0,010) 0,20) 0,04 | 0,058 | -0,11 | 0,08 | 0,0 | -0,2 | 7,1) 10,0 


37} 50 23) 0,014) 0,70) 0,004 | 0,24 | 0,049 | 0,063 | -0,06 |-0,02 |-0,1 | -0,2| 5,4) 7,1 


38 50 29) 0,009) 0,79) 0,007) 0,20) 0,04 | 0,055 | -0,16 | 0,04 |-0,2 | -0,2 | 6,3) 8,5 


39 50, 24) 0,013) 0,71  -0,023 | 0,23 | 0,049 | 0,058 | -0,07 | 0,10 | 0,21 0,1) 5,4) 6,4 


40} 50, 33] 0,008 | 0,80 | -0,008 | 0,21 | 0,041 | 0,051 | -0,07 | 0,17 | 0,0] 00, 6,6] 8,0 


41 50) 24) 0,013) 0,72 | -0,005 | 0,23 | 0,047 | 0,065 | -0,13 | 0,19 |-0,1 | -0,3| 5,6) 7,3 


42) 30, 21) 0,023) 0,63] 0,011) 0,19) 0,05 | 0,066 | -0,01 |-0,02 | 0,0) -0,3) 4,2) 5,4 


43) 30) 21) 0,023) 0,62 | 0,012) 0,18 | 0,052 | 0,069 | -0,04 |-0,02 | 0,1 | -0,3| 4,4) 5,8 


44) 30, 21) 0,023) 0,62 | 0,013) 0,19) 0,05) 0,07 | -0,09 | 0,01 0,1 | -0,4] 41) 5,7 


45] 30, 21) 0,023) 0,63 | -0,006| 0,17) 0,05 | 0,065 | -0,10 | 0,08 | 0,1 | -0,3, 43) 5,6 


46) 30, 32) 0,018) 0,66 | 0,014) 0,19 | 0,047 | 0,07 | -0,06 | 0,06, 0,1 | -0,2, 43) 6,6 


47) 30, 42) 0,016) 0,68 | 0,014) 0,18 | 0,042 | 0,064 | -0,06 | 0,00 | 0,2 | -0,2] 41) 6,2 


48 | 30) 46) 0,013) 0,72 | -0,006 | 0,17 | 0,045 | 0,061 | -0,17 | 0,11 0,1 | -0,2] 5,1) 6,8 


49 30) 51] 0,011) 0,75) 0,011) 0,16 | 0,041 | 0,065 | -0,18 | 0,18 | 0,2 | -0,1| 5,1) 7,8 


50 30; 31| 0,015) 0,69) 0,018) 0,13) 0,04 | 0,066 | -0,07 | 0,20) 0,3 | -0,2| 4,2) 6,6 


9.5 Empiryczne funkcje ortogonalne - POD 


Do analizy struktur przepływu zastosowano w pracy metodę rozkładu fluktuacji pola pręd- 
kości w bazie funkcji ortogonalnych (POD — Proper Orthogonal Decomposition). Analiza taka 
pozwala na zbadanie amplitud tych funkcji i wyznaczenie udziału w przepływie poszczegól- 


nych struktur turbulencji 


Metoda POD, podobnie jak na przykład analiza Fouriera, dokonuje aproksymacji danych 
przez ich rozkład w pewnej ortogonalnej bazie. Dla pewnego pola u(x,t) zdefiniowanego na 
pewnej podprzestrzeni (2 szuka się aproksymacji w postaci: 

K 


u(x,t)=Q,a'(r)g'(x), (9.31) 


k 
gdzie a‘ - współczynniki skalarne, a p*(x) wektory bazy spełniające warunek ortogonalno- 


Sci: 


J oh (x)? (x)d x=6, 4. (9.32) 


Q 
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W analizie Fouriera rozkładu dokonuje się w zadanej bazie funkcji trygonometrycznych. 
W omawianej tutaj metodzie empiryczne funkcje ortogonalne stanowią bazę, która jest 
wyznaczane na podstawie danych. Empiryczne funkcje ortogonalne są w badaniu turbulencji 
interpretowane jako struktury koherentne. Lumley ([45],[46]) zaproponował definicję struktu- 
ry koherentnej przez funkcje które zawierają maksimum energii. Są one liniowymi kombina- 


cjami o(x) , które maksymalizują następujące wyrażenie: 


arg max(|o(x),u(x,t)|), (9.33) 


o(x) 


gdzie wyrażenie (fg) oznacza tutaj iloczyn skalarny 


J JgdQ (9.34) 


a wyrażenie <> uśrednianie po czasie: 


T. 
(a(t))=lim + f a(e)at. (9.35) 
T>w T 0 


Jeśli o(x) maksymalizuje równanie 9.31 to oznacza, że pole rzutowane na o(x) ma 


największą energię w porównaniu do rzutowania na jakąkolwiek inną strukturę. 


W praktyce w metodzie POD analizuje się zbiór N pomiarów chwilowych u,„(x) fluktuują- 
cego pola prędkości. pomiary są wykonywane w różnych chwilach czasu. 
WAX) =a et (9.36) 
Zagadnienie maksymalizacji (równanie 9.31) można wtedy zapisać w postaci: 


\ N 


arg max = (0 (x),u,(x)) (9.37) 


a(x) n=l 


Rozwiązanie sprowadza się do zagadnienia własnego: 


RC=CA, (9.38) 


C - jest macierzą wektorów własnych, A - macierzą wartości własnych. 


R -jest macierzą kowariancji danych R=F" F , gdzie F jest macierzą danych: 
J J y g J 


u,(x,) U (Xe) u,(xp) 
p= u,(x,) uy (xy) A, lee (9.39) 
uy (xy) Uy (9) e uy (Xp) 


Macierz A zawiera wartości własne, które są ułożone na diagonali w kolejności rosnącej. 
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Kolumny macierzy C są wektorami własnymi — szukanymi funkcjami ortogonalnymi. 
Macierz ma następującą własność: 
C"CECC'2L, (9.40) 
gdzie I to identyczność. 
Wynika z tego, że otrzymane wektory są nieskorelowane, inaczej mówiąc są ortogonalne 


względem siebie, stąd nazwa techniki. 


Aby zbadać jak na przykład pierwszy mod ewoluuje w czasie, wystarczy obliczyć wyraże- 
nie: 
d,=FC,, (9.41) 
gdzie C, jest wektorem własnym odpowiadającym pierwszej wielkości własnej. Wektor d, 
jest wektorem o długości N odpowiadającej liczbie modów rozkładu. W ogólności możemy 


dla każdego Z, policzyć d, - główne składowe czasowe. 


Macierz danych może być zrekonstruowana przez obliczenie: 


P 
P=) 4 ¢,. (9.42) 


ial 


9.6 Wydruk programu PIV-Kor 


Program główny 
function [dx,dy]=PIV Kor(iml,im2,mask) 


% Program wykonuje algorytm PIV wraz ze znieksztalcaniem obrazow i 
filtrowaniem 

iml - pierwszy obraz 

im2 ="drugi obraz 

mask - maska, gdy chcemy niektore elementy obrazu wykluczyc z 
obliczen 


ale olo olo 


im1l=(double(iml)); 
im2= (double (im2)); 


[N,M]=size(iml); 


usuwamy ruch glowny 
[iml,im2,SX,SY]=odetnij (iml,im2); 


a= [65, 33, 25, 17, 15, 13]; % wektor zawierajacy rozmiary probki 
dla kazdej iteracji 
rmax= [25, 9, 5, 3, 2 , 1]; % wektor zawierajacy maksymalne 


dozwolone przemieszczenie dla kazdej iteracji 
odl= [31, 16, 10, sA 5 , 2]; % wektor zawierajacy rozdzielczosc 
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siatki dla kazdej iteracji 


media=[ 1, 1, 1, 1, 1, 1]; % wektor zawierajacy informaje o 
wykonaniu filtru dla kazdej iteracji 
wyg= [ 1, 1, 1, 0, 1, 0]; % wektor zawierajacy informacje o 


wykonaniu wygladzania dla kazdej iteracji 


[n,m]=size(iml); 
[ 
=X(ones(1,n),:); 


[ 
=Y(ones(1,m),:)'; 


@ Xw, Yw - polozenia punktow po deformacji 
% ustalenie wartosci poczatkowej: 

XW=X; 

Yw=Y; 

% imw2 - obraz drugi po deformacji 
im2w=im2; 


% petla po wszystkich wartosciach a 
for i=l:length (a) 


% wykonanie programu szukajacego przemieszczen i deformujacego obraz 
[Xw,Yw,im2żw]=PIVcor3(iml,im2,im2żw,Xw,Yw,mask,a(i),rmax(i),odl(i)); 

% obliczenie calkowitego przemieszczenia 

dx=Xw-X; 

dy=Yw-Y; 


wielokrotne filtrowanie 
if media (i) 


[xx, yy] =Siatka(iml,rmax(i),a(i) , odl(i)); 
XX=xx(1,:); 
yy=yy(:,1); 
dx=dx (yy, XX) ; 
dy=dy (yy, Xx) ; 


sx,sy]=median filter (dx, dy, ones (size (dx))); 
sx,sy]=median filter(sx,sy,ones(size(dx))); 
sx,sy]=median filter(sx,sy,ones(size(dx))); 
sx,sy]=median filter(sx,sy,ones(size(dx))); 
sx,sy]=median filter(sx,sy,ones(size(dx))); 
sx,sy]=median filter(sx,sy,ones(size(dx))); 
sx,sy]=median filter(sx,sy,ones(size(dx))); 
sx,sy]=median filter(sx,sy,ones(size(dx))); 


dx=interpoluj3(sx,xx,yy,m,n) .*mask; 
dy=interpoluj3(sy,xx,yy,m,n) .*mask; 


end 

% wygladzanie 

if wyg(i) 

rad=floor(1.7*odl(i)); 

rr=sqrt (repmat([-rad:rad],2*rad+1,1).*2+repmat ([- 
rad:rad]',1,2*rad+1).*2)<=rad; 

ss=sum(rr(:)); 


dx=kernl (dx,rr)/SSs; 
dy=kernl (dy,rr)/ss; 
end 


Xw=dx+X; 
Yw=dy+Y; 


end 


dodajemy przemieszczenia glowne odjete na poczatku 


dx=dx+SX*ones (size(dx)); 
dy=dy+SY*ones(size(dy)); 


spprzesuwamy dane, zeby zachowac ulozenie jak przed odcieciem 


if SX<0 
[n,m]=size (dx) 
dx=[zeros(n,SX),dx]; 
dy=[zeros (n,SX) ,dyl; 
end 


if SY<0 
[n,m]=size (dx) 
dx=[zeros(SY,m);dx]; 
dy=[zeros (SY,m) ;dy]; 
end 
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Program realizujący zniekształcenia obrazu 


function [Xw,Yw,Im2w]=PIVcor3(iml,im2,im2w,xw,yw,mask,a,rmax,odl) 


Program do PIV ze znieksztalceniami 
imi; im2 = oryginalne obrazy 


olo olo olo olo 


oo oo 


Im2w - wynikowy obraz znieksztalcony 


yw sa wspolrzednymi nieznieksztalconymi 


mask - macierz maski 

a wielkosc probki 

rmax - promien przeszukiwania 

odl - odleglosc miedzy punktami siatki 


ao ojo o oo 


tworzymy nieznieksztalcona siatke 
[n,m]=size(iml); 


X=[1:m]; 

X=X (ones(1l,n),:); 
Y=[1:n]; 

Y=Y (ones(1,m),:); 
Y=Y'; 

$$PIV 
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im2w - obraz znieksztalcony w poprzednim kroku 
xw, yw — znieksztalcone wspolrzedne z poprzedniego kroku 


XW,YW - wynikowe wspolrzedne znieksztalcone 


% jesli program jest wykonywany w pierwszym kroku to im2w=im2, a xw i 
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[XX, YY] =siatka(iml,rmax,a,odl); 
liczymy przemieszczenia 
[u,w,waga]=initmqd2 (im1, im2w, XX, YY,mask,a,rmax) ; 


SEeeee eee ee eye ee 


kis 


56666%% interpolacja na wszystkie piksele 
=interpoluj3(u,XX, YY,m,n) .*mask; 
(w, XX, YY,m,n) 


TT 
B 
3 
ct 
0 
R 
O 
O 
a 
u 
w 


QOQ2Q0200O0O0COOOOO0OOOOOC COO. 


Xw=interp2 (X,Y,xw,X+U,Y+W, 'linear'); 
Yw=interp2 (X,Y,yw,X+U,Y+W, 'linear'); 


% zgodnie z tym znieksztalcamy im2 otrzymujac imw2 


Im2w = griddata(Xw,Yw,im2,X,Y,'nearest'); 
Im2w = interpoluj_nearest (Xw,Yw,im2,X,Y,'nearest'); 
Im2w=interp2 (X,Y,im2,Xw,Yw, 'linear'); 


2 
Oo 
2 
6 


Szukanie lokalnych przemieszczeń 
function [dx,dy,w]=initmqd2 (iml,im2,X,Y,mask,a,rmax) 


stosowana do algorytmu ze znieksztalceniami 
anie podaje sie poprzedniego pola predkosci 
maska jest wpelni wymiarowa 


[n,m]=size(X); 

if ((n>1) & (m>1)) 
X=X(1,l:end); 
Y=Y(l:end,1); 

end 


m=length (X); 
n=length (Y); 


dx=zeros (n,m); 

dy=zeros (n,m); 

w=dx; 

sk2=sumowanie(im2.*°2,a); macierz sum kwadratow pikseli 
sk=sumowanie(im2,a); macierz sum pikseli 


émask=ones (n,m); 
for i=l:m @iteracje wzdluz x 


for j=l:n 5iteracje wzdluz y 
ssprawdzamy czy w danym polu w ogole cos jest 
if mask(Y(j),X(i))== 
continue 

else 

tworzymy powierzchnie korelacji 
kor=pow korf4(iml,im2,X(i),Y(j),sk2,sk,a,rmax); 
$kor=pow korf3(iml,im2,X(i),Y(j),a,rmax); 
[dx(j,i),dy(j,i),w(j,i)l=znajdz5 (kor); 


oo 
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end 
end 
end 


Tworzenie macierzy współczynnika korelacji 


function kor=pow korf4(iml,im2,x,y,s2k,sk,a,rmax) 


Tworzy powierzchnie korelacji wycinajac okno o boku a i srodku(x,y) z 
iml, nastepnie przesuwa je 
sna tle im2 maksymalnie o odleglosc rmax od punktu (x,y) 


wer=kwadrat (iml,a,x,y); %wyciecie probki 
wer=wer-mean(wer(:)); todjecie sredniej 


$poszerzenie probki do wymiarow okna przeszukiwan 
wls=zeros (at2*rmax) ; 

[s1,s2]=size(wer); 

wls(1:s1,1:s2)=wer; 


ale ole alo 


w2=kwadrat (im2,a+2*rmax,x,y); 


wls=zeros(size(w2)); 
[s1,s2]=size (wer); 
w1s(1:s1,1:s2)=wer; 


fl=fft2(wl1s); 
f2=fft2(w2); 


licznik 
kor=real((ifft2(£f2.*conj(£1)))); 


mianownik 
mil=sum(sum(wer.*2)); 


su2k=kwadrat (s2k,a+2*rmax,x,y); 
su2=kwadrat (sk, at+t2*rmax,x,y); 


mi2=(su2k - (su2.%°2)/(a%2)); 


mi=sqrt (mil*mi2) ; 


mi0=(mi==0); tsprawdzamy gdzie mianownik osiaga zero 
%laczenie 


kor=(kor./(mi0+mi)).*(-mi0); Suwazamy zeby nie bylo dzielenia przez zero 
kor=kor(1:2*rmax+1,l:2*rmax+1); 


Wyszukiwanie maksimum macierzy współczynnika korelacji 
function [d2,d1,maks]=znajdz5 (kor) 


% znajduje maksimum dla macierzy kor i dopasowuje do niego funkcje 
kwadratowa 


kor2=kor.*extrema (kor); 
maks=max (max (kor2)); 


[n,m]=size (kor); 
if maks== 
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Yes 

d1=0; 

d2=0; 
else 


[dl,d2]=find(kor2==maks) ; 
dl=round (mean(dl)); 
d2=round (mean (d2)); 


dl0=dkwadratowa (kor(dl-1,d2),kor(dl,d2),kor(dl+1,d2)); 
d20=dkwadratowa (kor(dl,d2-1),kor(dl,d2),kor(dl,d2+1)); 


dl=d1- (n-1)/2-1+d10; 
d2=d2- (m-1) /2-1+d20; 


end 
function ext=extrema (k) 


% dla pola k tworzy cos w rodzaju pola drugich pochodnych, przyblizajac 
je roznicami 
% sasiadujacych punktow 


[m,n]=size(k); 


a=k(2:m-1,2:n-1); 
b=ones (m-2,n-2); 
for i=-1:1 
fer j=-1i1 
if ( i==0 & j==0) 
else 
c=k(2+i:m-1+i,2+j:n-1+j); 
d=a-c; 
d=d>=0; 
b=b.*d; 
end 
end 
end 
ext=zeros (m,n); 
ext (2:m-1,2:n-1)=b; 


function d0=dkwadratowa (yl, y2,y3) 
m=2* (2*y2-yl-y3); 
if (yl==y3) | (m==0) 
H2"; 
d0=0; 
else 
d0=(y3-y1l) / (m); 
end 


Interpolacja 
function U=interpoluj3(u,X,Y,m,n) 
sinterpoluje pole predkosci zdefiniowane na danej siatce (X,Y)na 
powierzchnie 
calego obrazka o rozmiarach m x n 


U=zeros (n,m); 


X=int32 (X); 
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Y=int32 (Y); 


[gn, gm]=size (X); 

jezeli siatka zostala podana w formie macierzy to przerabiamy ja na 
wektory 

if ((gn>1) & (gm>1)) 

=X (1, 1:end); 

=Y (1:end, 1); 

end 


gna poczatek brzegi 


for j=1: (length (Y)-1) lewy 
us=itner(1,Y(j),X(1),Y(j+1),0,u(j,1),u(j+1,1),0); 
U(Y(J):Y(j+1),1:X(1))=us; 
end 


for i=1:(length(X)-1) *t@gorny 
us=itner(X(i),1,X(i+1),Y(1),0,0,u(1,i+1),u(1,i)); 
U(1:Y(1),X(i):X(i+1))=us; 

end 


for j=l:(length(Y)-1) prawy 
us=itner(X(end),Y(j),m,Y(j+1),u(j,end),0,0,u(j+1,end)); 
U(Y(J):Y(j+1),X(end) :m)=us; 

end 


for i=1:(length(X)-1) dolny 
us=itner(X(i),Y(end),X(i+1),n,u(end,i),u(end,i+1),0,0); 
U(Y(end):n,X(i):X(i+1))=us; 


end 
énarozniki 
U(1: sa 7 X(1))=itner(1,1,X(1),Y(1),0,0,u(1,1),0); 
U(1 do) ee SA BRERA 0,u(l,end)); 
U(Y poe X (end) :m) =itner (X(end),Y(end),m,n,u(end, ci +0; 0) 3 
U(Y (end) :n, ce X(1))=itner(1,Y(end),X(1),n,0,u(end,1) „0); 


srodek 
for i=l:length(X)-1 
for j=l:length(Y) - 
us=itner(X(i a (j),X(i+1) ,Y(j+1) ,u(j,i),u(j,i+1) ,u(j+1,i+1) ,u(j 


+1,1)); 


U(Y(J):Y(J+1),X(i):X(i+1) ) =us; 


end 


Filtr medianowy 
function [sx,sy]=median filter(vx,vy,mask) 
% filtr medianowy dla pola wektorowego 


[n,m]=size(vx); 
emask=mask (Y,X); 


x=repmat ([1:m],n,l); 
y=repmat ([1:n]',1,m); 
x=x (logical (mask)); 
y=y (logical (mask) ); 
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vx=vx (logical (mask)); 
vy=vy (logical (mask)); 


lvec=8; %liczba wektorow 


sx=zeros (n,m); 
sy=zeros (n,m); 


for i=1:length (x) 
% szukamy 8 najblizszych sasiadow 
[R,I] =sort((x-x(i)).*2+(y-y(1)).72); 
Vx=vx (I (1:lvec+l)); 
Vy=vy (I (1:lvec+1)); 


@ tworzymy z wektorow macierz kwadratowa 
dVx=repmat (Vx',lvec+l,1); 
dVy=repmat(Vy',lvec+l,l1); 


żobliczamy roznice 

D=sqrt ( (dVx-dVx') .*%2+(dVy-dVy') .*2); 
md=sum (D) /lvec; 

med=find (md==min (md) ) ; 

med=med (1); 


delta med=md (1,med) ; 


if D(1,med)>delta med 


sx(y(i),x(i))=Vx (med); 

sy(y(i),x(i))=Vy (med); 
else 

sx(y(i),x(i))=vx(i); 

sy(y(i),x(i))=vy(i); 


end 


end 


Tworzenie siatki punktow pomiarowych 
function [X,Y]=siatka(Im,rmax,a , odl) 
% Tworzy siatke dla obrazka Im uwzgledniajac maksymalne przesuniecie 


rmax i szerokosc okna, 
% w odleglosci odl 


X=[((a-1)/2+rmax+1) :odl: (n-(a-1)/2-rmax)]; 
Y=[((a-1)/2+rmax+1) :odl: (m-(a-1)/2-rmax)]; 


X=repmat (X,n,1); 
Y=repmat (Y',1,m); 
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Korelacja obrazu z zadanym wzorem przy użyciu FFT 
function a=kernl (im, kern) 


im - macierz wejsciowa 
kern - wzor - musi byc mniejszy od im 


ojo Alo 


tworzymy macierz o takim samym wymiarze jak przetwarzana macierz im 
zawierajaca wzor kern 
reszte wyrazow zerujemy 

n,m]=size(kern); 

mask=zeros(size(im)); 

mask(1:n,1:m)=kern; 


™ go o olo 


3 liczymy korelacje 
a=real (ifft2 (fft2 (double (im) ) .*conj(fft2 (mask)))); 


wyrzucamy niepotrzebny kawalek 
a=a(l:end-n+l,l:end-m+1); 


A=zeros(size(im)); 
A(floor(n/2)+1l:end-floor(n/2),floor(m/2)+1:end-floor (m/2))=a; 
a=A; 


Wycinanie kwadratowego fragmentu z obrazu 
function im=kwadrat (Im,a,x,y) 


Im - macierz wejsciowa 
a = dlugosc boku kwadratu 
x,y wspolrzedne 


alo olo olo 


[n,m]=size (Im); 


% funkcja window wycina z macierzy Im kwadrat o boku a i srodku x,y 


if floor(c)~=c 
t Blad = a musi byc nieparzyste"! 


else 
starty=(y-c)*((y-c)>=1)+((y-c)<1); 
startx=(x-c) * ((x-c) >=1)+((x-c)<1); 
endy= (y+c) * ( (y+c) <=n) +n* ((y+c) >n) ; 
endx= (xtc) * ((x+c) <=m) +m* ( (x+c) >m) ; 


im=Im( starty:endy , startx:endx ); 
end 
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Lokalne sumowanie pikseli po zadanym obszarze 
function s=sumowanie (im,a) 


program liczy dla kazdego punktu obrazu sume pikseli z jego otoczenia 
znajdujacych sie w kwadracie o boku a 
im - macierz wejsciowa 
a - szerokosc kwadratu 


ale alo olo olo 


[n,m]=size(im); 
mask=zeros (n,m); 
mask(l:a,l:a)=l; 


d=real (ifft2((fft2(im)).*conj (f£t2 (mask) ))); 
s=d; 
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